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PROJECT SUMMARY 


PURPOSE: 

To improve the quality and reliability of numerical simulations in 
computational mechanics. 


AREAS OF RESEARCH: 

Error estimation, solution enhancement, and mesh conditioning for 
finite element, finite volume, and finite difference meshes and 
solutions for problems in linear elasticity, steady state heat 
conduction, viscous incompressible flow, viscous and inviscid 
compressible flow. 


RESEARCH RESULT: 

An operational software package, AUDITOR, incorporating the 
following advanced features: 

• Accepts unstructured 2D and 3D neutral files for finite 
element, finite volume, and finite difference meshes of arbitrary 
spectral order; 

• Incorporates an object-based data structure; 

• Pre-processes meshes for conditioning characteristics; 

• Computes local and global solution error with advanced 
residual error estimation techniques; 

• Post-processes initial solutions to derive enhanced, higher 
order results; 

• Provides advanced 3D visualization; 

• Operates in a modern, X-MOTIF graphical user interface. 

POTENTIAL APPLICATIONS: 

The AUDITOR code may be operated as a standalone adjunct to 
many commonly used simulation codes, or its features may be 
integrated into such codes, thereby enhancing their reliability 
and robustness. 
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1 Introduction 


This report is the final report for the Phase II project entitled ”Pre- and Postprocessing 
Techniques for Determining Goodness of Computational Meshes”. This project has focused 
on the development of a number of new ideas and methodologies which could significantly 
improve the quality and reliability of numerical simulations in computational mechanics. 
The methodologies are based on the concepts of error estimation, solution enhancement, and 
mesh optimization each of which plays a integral part in the success of the modeling effort 
and reliability of the numerical results. The ultimate goal of this effort is to provide a stand- 
alone pre- and postprocessing package which is capable of accepting arbitrary unstructured 
meshes and numerical solutions from various class of numerical solver techniques and provide 
graphical and tabulated feedback on the quality of the mesh and reliability of the solution. 

With this rather general set of objectives in mind, a Phase I research and development 
effort was undertaken aimed at both a-priori and a-posteriori improvement of the quality of 
numerical solutions for a large class of linear and nonlinear computational mechanics prob- 
lems. In particular, this effort focused on error estimation, solution enhancement, and mesh 
conditioning for finite element, finite volume, and finite difference computational meshes and 
solutions for linear elasticity, steady state heat conduction, viscous incompressible flow, and 
inviscid and viscous compressible flows. The specific objectives of the effort were as follows. 

• The development of two-dimensional and three-dimensional preprocessing techniques 
which provide the user with an a-priori indication of the quality of the initial compu- 
tational mesh. 

• The development of two-dimensional and three-dimensional a-posteriori error estima- 
tion techniques for accessing the cell-wise errors in user supplied numerical solutions. 
These estimates are based on residual error estimation techniques which provide both 
a local and global estimate of the error in the numerical solution. 

• Post-processing methods were to be developed for two-dimensional and three- 
dimensional simulations which provide for higher order solution extraction. 

• Graphics and visualization capabilities were to be developed to provide the user with 
a graphical interface to view the computational domain, the solution components, 
extracted enhanced solution data, and distributions of the solution error. 

This report contains detailed descriptions of the work done toward meeting all of these 
objectives during the course of the Phase II effort. In particular, a summary of the research 
and development progress is provided in Section 2. This is followed in Section 3 by an 
overview of the theory, algorithms, and methodologies studied in the course of this Phase II 
effort, many of which have been incorporated into the AUDITOR pie- and postprocessing 
package. In Section 4, an overview of the functionalities and features which comprise the 
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A UDITOR package are presented. This is followed in Section 5 by a summary of the results 
of several benchmark problems that have been solved to validate various components of 
the code. In particular, representative results for two-dimensional and three-dimensional 
heat conduction, linear elasticity, compressible and incompressible flows are presented here. 
Finally, Section 6 concludes this report with some observations about potential areas for 
continued research and development. 
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2 Summary of the Phase II Effort 

The principal objectives of the Phase II research and development effort have been outlined 
in Section 1. The ultimate goal of this project, however, has been to synthesize these 
features and functionalities into a fully operational tool for the two- and three-dimensional 
pre- and post-postprocessing of user supplied computational meshes and numerical solutions. 
The culmination of this effort has resulted in what we believe is a unique software package 
that contains features and analysis capabilities that are collectively not available in any 
other commercial or private analysis package. In particular, we believe that the residual 
error estimation capabilities, the ability to accept unstructured meshes and solutions from 
finite element, finite volume, and finite difference solvers, the scope of the class of problems 
accepted (linear elasticity, steady state heat conduction, compressible and incompressible 
fluid dynamics), and the ability to audit hp-adapted meshes, represents a significant advance 
in the area of computational mechanics. 

The technical accomplishments that have been reached during the Phase II effort are 
summarized as follows; 

• The development of residual error estimation capabilities which incorporate the fol- 
lowing features or functionalities 

— are applicable to finite element, finite volume, or finite difference numerical solu- 
tions 

— ^Aral /rpll.wisp^ trlobal estimates of the solution error in the energy 
norm 

— provide accurate error estimates of solution error for unstructured meshes which 
may be composed of cells with arbitrary spectral order 

— do not require the specification of boundary conditions; the appropriate boundary 
conditions are assumed to be satisfied by the approximate solutions. 

— solve a sequence of local problems to obtain estimates of the error and thus avoids 
the assembly and solution of a large system of equations 

— provide an interface to the graphics package to view the results of the error esti- 
mation process 

• The development of mesh postprocessing capabilities which incorporate the following 
features: 

— are applicable to finite element, finite volume, or finite difference numerical solu- 
tions 

— extract enhanced higher order results for derived quantities based on the initial 
solution 
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— provide and interface to the graphics package to view the results of the enhance- 
ment process 

• The development of mesh preprocessing capabilities for performing an initial mesh 
preconditioning analysis. In particular, the efforts in this area have resulted in the 
following: 

— a mesh preconditioning module which checks the smoothness and orthogonality 
of the initial mesh 

— a module which allows the user to interactively specify criteria on which to test 
the mesh characteristics 

— a diagnostic feed back both in a tabular form and graphically though the graph- 
ics/visualization package. 

• The development of motif graphical user interface which is mouse driven and is com- 
patible with X.window platform. 

• The development of various components of a three-dimensional graphics package for 
viewing unstructured meshes and solutions. Among the features which comprise this 
package are; 

— isosurface plotting for both volumes and surfaces 

— planar slicing 

— velocity vectors 

— xy-plotting capabilities 

— local error estimate plotting capabilities 

— mouse driven zoom, pan, rotate, and tumble features 

— postscript hard copy capabilities 

• The development of gridfile interface with the following options, 

— reads a gridfile in a predefined COMCO format (see the AUDITOR user manual 
for details) 

— reads a neutral file generated from P3/CFD or PHLEX 

In addition to these milestones significant progress has also been made in the following areas: 

• A generic neutral file reader has been developed which accepts unstructured two and 
three-dimensional computational meshes in a given format and directly translates this 
data into cell wise data that is loaded into the AUDITOR database. The scope of this 
package includes: 
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— a capability for accepting finite element, finite volume, or finite difference meshes 
and solutions in various formats 

— a capablility for accepting various types of solution data including solution deriva- 
tives at arbitrary specified locations 

— a capability for converting various element types and associated data formats 
• The development of an extended object based data structure which includes: 

— dynamic memory management 

— objects structures and methods 

— stacks, heaps, and queues 

— manipulates multiple element types 

* hexahedral elements 

* tetrahedral elements 

* prism elements 

* shell elements 

The final code incorporates some features that are common to other software packages 
that are under development at COMCO. These include the current data structure on which 
AUDITOR is built, and the graphics and postprocessing capabilities. In addition, initial 
support for Phase III productization of some of this technology in the area of computational 
fluid dynamics has been received by a partner company and couid resuit in a commercial 
release in a new residual error estimation feature for P3/CFD within the next six to twelve 
months. 
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3 Technical Approach 


This section describes the technical aspects of the methods used in A UDITOR to perform the 
various tasks. Also included in this section are summaries of the research studies performed 
during the project in developing these methods. 

3.1 Mesh Preprocessing 

The mesh conditioning module analyzes the mesh, quantifying its various attributes and 
searching for inconsistencies. This analysis is usually done by evaluating special functionals 
designed to quantify various desirable features of a mesh, such as smoothness and orthogo- 
nality. The definition of these functionals varies from author to author, but all of them are 
invariably dependent upon the geometric features of the mesh. 

For example, Carcaillet, et. al. [1] use functionals based upon vectors oriented along cell 
edges. Smoothness is defined according to 

#nodes #edges(i ) 

F. Mh = £ £ ||r„|| ! 

t=l J = 1 

where i is the set of corner nodes, j is the set of edges connected to node i and ||r tJ || is 
the length of a particular edge. Orthogonality is measured according to the functional 

# nodes #faces(i) 

^orth = y! y (**01 ' ** 02 ) 

*=1 i=i 

where i is the set of corner nodes, j is the set of faces connected to node i and r,ji and 
r tJ 2 are the two edges on face j connected to node i. 

Other authors, such as Jacquotte [2,3], have defined functionals based on cell deformation 
ideas. Specifically, these functionals use the invariants of the left Cauchy-Green tensor of the 
deformation of a reference cell into an actual cell. In two-dimensions, this functional takes 
on the form 


F = C 1 { /, - 2J) + C 2 (J - l) 2 

and in three-dimensions 


F = C 1 (h + 1 2 -6J) + C 2 (J-l) 2 

where 7 X and / 2 are the invariants of the left Cauchy-Green tensor C of the transformation 
from a reference cell into a physical cell x = *(£). J is the Jacobian of the transformation 
and Ci and C 2 are constants used to weight the two terms. These quantities can be written 
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I\ — tr C 
I 2 = trCofC 
J = det Vx 

In these functionals, the first term measures the “conformality” of the deformation while the 
second term is a volume control term. Together, these two terms measure the smoothness 
and orthogonality of the mesh. 

The most common technique for measuring smoothness and orthogonality, however, was 
introduced by Brackbill and Saltzman [4]. This method uses the transformation from a 
reference (£, 77 ) to physical space (x,y) to define functionals for smoothness 


■^smooth £[(V£-V0 + (Vi 7 . Vy)]d(dy 

Jn 

and orthogonality 

h = /.l(V£ • 

Jn 

Most authors, such as Kennon [5] and Demkowicz [ 6 ] however, write these functionals in 
terms of the transformation * = »(£) 


smooth 


Ja 


(dx*\ . (dy*\ . (dx*\ . (dy*\] 1 

[ysa r w r\* r m n 


and 


F orth 


dxdy dx dy ^ 2 


d£ dy dy J J 


d(dy 


As the work of these authors indicates, the condition of the mesh can be measured by 
examining the transformations from the reference cell to the physical cells. In particular, this 
can be done by evaluating the jacobians of the transformation from the master to physical 
element for each element in the mesh. The values of the jacobians and their determinants 
can be used to indicate several different measures of mesh quality. For example, a negative 
determinant indicates an element which is “inside out” or at least partially folded over on 
itself. Wide variations in the values of the determinant in an element imply a highly distorted 
element with very small and/or large angles which often result in a poor local approximation. 
Furthermore, a non-constant determinant usually results in a certain amount of error when 
quadrature rules are used to integrate over the element. This error may become significant 
the more the determinant deviates from a constant especially if it becomes very small. 
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As a result, mesh analysis based upon the element jacobians would identify highly dis- 
torted elements and elements with bad (negative) determinants, and then allow the user to 
view the offending elements. Furthermore, the analysis would evaluate statistical measures 
of the variations of the jacobians in order to provide an overall indicator of the mesh quality. 


3.2 Error Estimation Techniques and Research 

A significant portion of the work performed on this project was spent researching error 
estimation techniques which would be suitable for the wide range of problems under con- 
sideration (linear convective heat transfer through compressible viscous flow). This section 
summarizes these efforts following a brief introduction on the purposes of error estimation. 

3.2.1 Why Bother to Estimate Errors? 

The first, foremost, and obvious reason an analyst would want to estimate errors in an ap- 
proximate solution is reliability. Simply put, if engineering design or analysis decisions are 
to be made based upon data produced by an approximate solution technique, this technique 
must have some way of providing some assurance that the calculated data is accurate. Un- 
fortunately, most engineering analyses rely primarily upon the expertise of the user. While 
there is no substitute for experience, a user simply does not have enough information to 
make an informed judgment. Traditionally, when examining a solution, a user will look for 
inconsistencies in the solution or changes in the solution from one mesh to the next. For 
example, in an elasticity problem, a user might look for discontinuities in the computed 
stresses or changes in the strain energy of the computed solution from one mesh to the next. 
(Some analysis programs refer to this change in the energy of the solution as the “error in 
energy.”) While these techniques indicate when error is present, the converse is not true. A 
computed solution may have continuous derivatives and may not change significantly from 
mesh to mesh but still be inaccurate. Consequently, these techniques used to analyze a com- 
puted solution have no direct itlationship with the actual error, they merely are features that 
the error may or may not exhibit. 

In reality, when attempting to estimate a solution’s reliability, the user wishes to know 
the difference between the computed solution and the solution which satisfies the governing 
differential equation and boundary conditions exactly. Furthermore, the user may not be 
interested in the solution but in some combination of its derivatives (such as stresses, heat 
fluxes, or skin friction forces) and may be interested in the values of these quantities only in 
certain places. Suppose the analytic solution were known. Figure 3.1 shows the computed 
solution, the error in the computed solution, and the errors in the solution’s x- and y- 
derivatives for a mesh of 18 quadratic elements where the analytic solution is a simple 
bi-cubic polynomial. A reliability analysis would want to know how good the solution is 
in some region of the mesh. As Fig. 3.1 clearly demonstrates, even with knowledge of the 
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analytic error, a reliability analysis can still be difficult to perform. 

For this reason and others, error estimates are generally calculated in terms of a norm of 
the estimated error which in turn is usually “normalized” by the same norm of the solution. 
For example, if the L 2 norm is being used, the error norm for a particular computational 
cell, fl* would be: 



where u and it/, represent the analytic and computed solutions, respectively. In this form, 
the error norm represents an “average” error over the computational cell as a percentage of 
the solution (i.e., an L 2 error norm of 0.05 would imply that the computed solution is on the 
average about 5 percent different from the analytic solution). 

In a typical error estimation analysis, the above error function, u - u h , is replaced by 
an estimate, e/,, and the estimate is normalized by the norm of the computed solution. 
Additionally, since the user is often interested in the derivatives of the computed solution, 
the Hi norm or the “energy” norm is used instead of the L 2 norm. 

Notice that in order to compute an error estimate norm, the estimated error must first 
be calculated. Under certain circumstances, this function estimating the error can be added 
to the computed solution resulting in a new solution with greater accuracy. 

In addition to determining the reliability of a computed solution, error estimates provide 
information to the user which can be used to change the mesh to improve the solution. This 
idea is based upon the premise that replacing a cell that has a large error with a number of 
smaller cells (or a cell with a higher spectral order) will reduce the error in that region of the 
mesh. For diffusion dominated problems, this premise is rarely violated. For problems with 
significant convective phenomena, however, the errors in elements which are “upstream of 
the element must also be considered. 


3.2.2 The Dual Norm and Residual Error Estimates 

Perhaps the most natural functional setting for the finite element method is to assume that 
solution space U is a subspace of if 1 , whereas space V is a topological dual of U, V = U' . 
This requires that the residual is to be measured using the dual space norm 


INI u< = su p 

v£U 

v?0 


IMI 


( 1 ) 


This idea has led to the element residual method which will be presented in this and the 
following sections. Furthermore, this technique was carefully examined in a number of papers 
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PROJECT: Enor Demo 


FlNflE EIJiMENT SOLUTION 


PROJECT. Eiiof Dciiuj 


ERROR IN SOLUTION 



MIN- 0 131736 
MAX-0 1317361 


MIN— 0 023149 
MAX -0 0231417 



PROJECT; Enor Demo ERROR IN X DERIVATIVE 



PROJECT: Enor Demo ERROR IN Y DERIVATIVE 



Figure 3.1: Computed solution, error in computed solution, and error in solution’s derivatives 
for a mesh of 18 quadratic elements and a bi-cubic analytic solution. 
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(see [7,8,9]). 

Obviously the norm denoted by equation (1) depends upon the choice of the norm on 
the original space U. Considering for instance the linearized, steady-state Navier-Stokes 
equations 


— i/A« + (U • V)u + gradp = / 

div u = 0 


( 2 ) 


with (for simplicity) Dirichlet boundary conditions on 


u, we define the space U as 


U = {(u,p) (E Hl(Q) x L 2 (Sl)} 

and define the residuals corresponding to the two equations in (2) as 


(3) 


sup 

V€Hl({l),V*0 


sup 

qeL 2 (Q),q? 0 


X 


(i/'VuhVv + ( U ■ V)uhV — phdivv — fv)dx 

M 

/ divuhqdx 

H 


(4) 


With the choice of L 2 -norm for q, the second of the residuals residuals reduces just to the 
jf, 2 -residua.l of the divergence 


\\divu h \\ L i (n) (5) 

whereas the choice of the norm for v is neither obvious or unique. 

Once a norm is selected for v we may proceed by introducing local Neumann problems 
over elements f Ik as 


-vA<p K + {U -V)<p K = r h 

t(<Pk) = <*K9h 


( 6 ) 


where r/, is the element residual defined as 

*h = f + vAu h -(U • V)u k - gradp h (7) 

and g h is the flux residual including the jump in both the pressure (if it is discontinuous) and 
viscous fluxes (due to the discontinuity of derivatives of u h). The “distribution parameter” 
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ax accounts for distributing the flux residual to both neighboring elements in such a way 
that the local Neumann problems are well-posed. 

Once the local problems are solved the element error indicators are defined as 

Vk = \\<Pk\\ 2 + \\ divu h\\ 2 ( 8 ) 

This section focuses on the two fundamental issues: 

1. How to select the norm. 

2. How to determine the distribution parameters ax- 

We emphasize that the use of local Neumann problems does not require any orthogonality 
conditions for the approximate solution u ^ and therefore may be applied to resulting from 
virtually any approximate method. 

Error Estimation for Stokes’ Flow 

One central problem of the project was the development of a general “black box” error 
estimation package for assessing the accuracy of approximations obtained using finite ele- 
ment, finite difference, or finite volume techniques. Recent theoretical work [10] shows that 
this goal is now attainable and can be justified on a rigorous mathematical basis, albeit for 
certain classes of problem. 

Unfortunately, as yet, theory is not well developed for the estimation of errors in approx- 
imations to the incompressible Navier-Stokes equations: 


— vAu + (u • V)u -f Vp = f 
V • u = 0 


( 9 ) 


Ultimately, we shall produce an error estimation procedure for these equations. However, 
rather than attempting to deal with them in their entirety, we propose to focus our attention 
on the key aspects of this problem by examining several intermediate model problems. The 
key issues are: 


• In what norm should the error be measured? In particular it is necessary to decide on 
how to “balance” the error in each of the equations or variables u,p appearing in (9). 

• How should one deal with the nonlinear term u ■ Vu? The question here is not purely 
related to the nonlinearity but also to the nonsymmetric nature of the term. Once 
again the question of norms arises. In order to obtain a measure of the error which 
will be suitable for mesh refinement and enrichment purposes, it is necessary to weight 
the error upstream since refinements do not only lead to improvement locally but also 
have an impact downstream. 
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With these considerations in mind, we shall focus on the following two model problems: 
Stokes ’ Flow 


— i/Au -f Vp = / 

► 

V u = 0 

and Heat Flow With Advection 


( 10 ) 


— i/Ait + (a • V)u = f 


( 11 ) 


The purpose of studying (10) is to investigate ways of “balancing” errors between the two 
entities, and by studying (11) one may gain insight into what norm can be used to weight 
upstream errors correctly. 

Stokes’ Flow 

Let us initially suppose that Uh and ph are some approximate solution of the Stokes’ 
Problem: 



—An + Vp = f } 

i 

V • u = 0 J 

y in J7 

(12) 

subject to 

u = 0 on dtl 

(13) 

and 

J n pd X = 0 


(14) 

We denote 





A' = Hi(Sl) x Hi(Sl) 

M = £2(S1) = {« € i ! ((l) : Jqdx=o\ 


( 15 ) 


The variational form of (12) then becomes: 
Find (tt,p) € X x M such that 
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( 16 ) 


where 


B(u,v) + b(v,p) = (/,t 0 V v 6 X 
b(u,q ) = 0 V q G M 


= J grad u: grad v dx j 
6(r,p) = -L pdivv dx 

(/,») = jjv 


1 dx 


As a preliminary assumption we assume 


w/i € X , p* 6 M 

but that otherwise and pk are arbitrary. 

Denote e = u — Uh E X and E = p Ph € M • Then fiom (16) we 
(e,E) € X x M : 


£(e,t>) + 6(v,£) = (/,»)-%,t))-H»,Pk)V®6^ 1 
6(e,g) = -b(u h ,q) V q e M J 

Balancing of Errors 

Let e > 0 be some user-specified parameter, the^ interpretation of which 
apparent later in this discussion. Then we define (e, Ej € X x M: 


B(e,v) + b(v,E ) = {f,v) - B(u h ,v) - b(v,p h ) V v e X 

-e(E,q) + b(e,q) = -b(u h ,q) VqtM 

Since v € X then V ■ v e M and so choosing q = V • v in (20) we obtain 

- e(E,V-v) + b(e,V-v) = - b(u h ,V-v ) or 
eb(v,E ) - (V • e, V • r) = (V-v h ,V-v) 


(17) 


(18) 

obtain that 


(19) 

will become 


( 20 ) 
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( 21 ) 



and now substituting (21) in (20) gives 


e € X: B(e,v) + ~{V-e,V-v) = (/,v) - B(u h ,v) - b(v,p h ) 


-(V.u fc ,V.») 


The bilinear form appearing in (22) may be written as 


a(e,v) = B (e,v) + - (V • e, V • v) 


where 


= [ (Se) t D(Sv)dx 

Jo 


i • 


5 = 0 


0 i. 

oy 

d_ d_ 

dy dx 


1 -(■ £ 1 — £ 0 

D = - 1-el+eO (25) 

£ 

.0 0 £ . 

With this notation, (22) becomes 

c 6 X : a(e,v) = (f,v) - a(u h ,v) - b(v,p h ) V v € X (26) 

Formally, the problem (26) is a linear elasticity problem. That is, we can in principle solve 
(26) to obtain e « e. The natural norm for measuring the solution of (26) is 


||e||| = a (e, e) = B (e, e) + - || V • e||^ (n) (27) 

and it is with respect to this norm that we propose to assess the accuracy of w/, and indirectly 
p h . Notice that the user-specified parameter £ acts as a weight on the divergence of the error. 
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Thus by selecting an e, the user specifies how great an effect the violation of incompressibility 
has on the error. Thus, in order to deal with the Stokes problem it is necessary to first deal 
with the linear elasticity problem, which is of independent interest also. 

Linear Elasticity 

Let us consider the following problem: 

Find u such that 


— = f in ft, (T = DSu (28) 

subject to 

H<t = g on T n 

where we suppose the data g and f are compatible and 

• D is a general elasticity matrix 

• H is the matrix formed using the components of the unit outward pointing normal 
n = (n r , n y ) T as follows 


(29) 


H = 


Tig 0 Uy 

0 n y n x 


Let 


(30) 


X = H\tt)xH'{a) 

then the weak form of (28)-(29) is: 

Find u € X such that 


J (SuYD(Sv)dx = j f l vdx-\- J g l vds 


(31) 


for all v 6 X. Let Uh £ X denote some approximation to the elasticity problem (31) and 
denote e = u — u h . Then 


a(e,v) = (/,v) + (</,v) rw -a(w fc ,i;) V«GX (32) 

The problem in which we are interested is the estimation, and preferably bounding, of 
j||e||| = yja{eYe). In order to facilitate our analysis it is convenient to introduce a zeroth 
order term; for A > 0 define 
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( 33 ) 


a(u, v) = a(u,v) + (Au,v) 

The significance of A will be discussed later. We then let e € X: 


a (e,v) = a(e,v) V v € X 


(34) 


and consider the problem of estimating a (e, e). Obviously as A is chosen smaller and smaller, 
we find 


lim a (e, e) = a(e, e) (35) 

A-.0+ 

Localization 

In order to minimize the cost of obtaining the error estimators, we first of all break (32) 
into a sequence of local problems on each element. 

Let us define various quantities associated with the local problem on Oa' 

Ot (*) = fi x ) + S t DSuh (36) 

and 


fa — HDSuh on 50 a - 0 IV 

IMS) = < 37 > 

l -oklIH DSu h ] on rXt K n iJSl L 

where [•] denotes the jump in tractions of the approximate solution across the interelement 
boundary. Here a KL is a 3 x 3 diagonal matrix weighting of the approximate tractions from 
elements Ha' and Ha on the edge of Ta a = H BQ,l- The weighting should ideally be 
chosen so that 


(. HDSu h ) a : = a KL (. HDSu h ) K + a LK ( HDSu h ) L w H DSu (38) 


(i.e., the weighted average of the tractions) represents a good approximation to the true 
traction. Obviously, in order to obtain consistency we require 

OCKL + a LK = I ( 39 ) 
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where / is the 3 x 3 identity matrix. In most applications to date, the choice oikl — was 
sacrosanct, but it can be shown [10] that this choice has severe repercussions regarding the 
accuracy of the resulting estimator. The local problems then become 

Find <p K € H 1 (Qk) X H l {n K ) : for all w £ H l (Sl K ) x 

[ (SvKYDiS^dx = [ r K {xyu>dx + <f R K {sYu>(s)ds 

W j *k Jd ° K ( 40 ) 

- Sk I “ dx 

JQk 

where 

s « = \<h\ iL TK{x)dx + L k RK{s)ds } (41) 

The physical interpretation of the term involving 6k, which arises naturally from the math- 
ematical analysis, is to enforce equilibrium on the local problem. Having obtained <p K we 
then define 


vl=[ (Sv> K )'0(S*> A .)<te + 1^1 l«ftf ( 42 > 

The following result, generalizing a similar result in [10], can be shown: 


ll|S||| ! + A||e|| I <E’lt (43) 

Ok 

Of course this gives an upper bound on a (e, e). However, let us suppose that we can choose 
clkl (see (37)— (39)) in such a way that 6 K = 0. Then letting A -» 0 + and using (33), (42), 
and (43) gives 


ll|e|ir<E/ {S<PkY D(S<PK) dx (44) 

o K Jn « 

This result is perhaps rather astonishing, proclaiming an upper bound on the discretization 
error. The key is the assumption that S K can be driven to zero. This may be understood 
intuitively as follows. In order to estimate the error we can essentially only base our obser- 
vations on the residuals r K and R h However, there is a third type of residual, hitherto 
neglected namely that the true solution is in equilibrium locally. This corresponds to 6 K = 0, 
therefore we should incorporate 6 K into our analysis of the error. The above discussion shows 
how to accomplish this. 

Flux Splitting 
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We now make some remarks on how okl should be chosen to ensure 6 k = 0. The most 
crude, but simplest, way to achieve this is by solving the quadratic programming problem 

min^ |5 a-| 2 subject to cukl + ola = / on 50 a H 50 l (45) 

Ok 

This represents a computation on a global scale and is therefore rather expensive. However, 
a sensible strategy would be to choose okl = during the adaptive analysis and then 
perform (45) only once in the final phase when the final solution has been obtained. The main 
drawback of this technique is the expense. However, if we assume that the approximation 
arises from an h-p element computation then considerable savings result. 

Flux Splitting for h-p Finite Element Computations 

Let us now assume Uh, is an h-p finite element approximation. In [10] it was shown that 
oal for which 6k vanishes can be calculated by local computations only. This makes the 
cost very inexpensive and is suitable for parallelization, reducing the time to a negligible 
amount. We shall not give the precise details at this point but it is worth noting that the 
only computations required are the solutions of small (at most 6 x 6 in two dimensions) 
systems of algebraic equations. 

3.2.3 Error Estimation Using the Flux Balancing Technique 

In this section we provide a more complete discussion for the flux balancing technique for 
equilibrating the local problems during error estimation. For additional details on this tech- 
nique, we refer to Ainsworth and Oden [11]. Although the following discussion is presented 
in the context of finite element techniques, the approximate solution whose error is being 
estimated could have been calculated using almost any approximation technique. In this 
case, the finite elements merely represent computational cells. This flux balancing technique 
was implemented into a two-dimensional finite element test code and used to compute er- 
ror estimates on several classes of problems. After the following discussion, we present the 
results of sample computations performed for a vector-valued model problem. 

Self— Equilibrating Flux Splittings 

First we begin with some historical remarks from which the idea of self-equilibrating flux 
splittings evolved. Consider for a moment the Neumann problem: 

Find w such that 


dw 

— A w(x) = f(x ) in fi subject to -r— = g on 5fl 

on 


The well-known condition for the existence of solutions to this problem is that 

f f(x)dx + <f g(s)ds = 0 
Jn JdCl 


(46) 


( 47 ) 
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When considering the problem of obtaining a posteriori estimates of the error in a finite 
element approximation of (46), one is led to consider local Neumann problems on each 
element of the partition with data related to the residual / + A u and to the jumps in the 
normal fluxes across interelement boundaries. Of course it is necessary that the data satisfy 
the compatibility condition (47) in order for the local problem to be well posed. 

In general the compatibility condition will be violated. Various ad hoc ways of reformu- 
lating the local error problem have been used to circumvent this incompatibility problem. 
However, Kelly [12] proposed that the fluxes be split unequally between elements in such a 
way that the compatibility condition be satisfied, referring to this idea as self-equilibration. 

In one dimension for a linear approximation, the splitting could be found explicitly in 
terms of the interelement jumps in the normal fluxes and the residual integrated against the 
basis functions. However, in higher dimensions such an explicit formula could not be found 
and instead Kelly resorted to minimizing the sum over all the elements of the squares of the 
amounts by which the compatibility condition is violated. It was found that splittings could 
be computed for which the lack of equilibrium was negligible. 

Kelly [12] gives numerical results demonstrating the significant improvement in perfor- 
mance of the error estimators when the self-equilibration is applied. However, the idea has 
not been widely adopted by the finite element community, in spite of the increased emphasis 
on a posteriori error estimation. The primary reasons are probably because the procedure 
is extremely costly (due to a global minimization) and is needed only for operators such as 
the Laplacian. 

Generalized Flux Splitting 

More recently, once again working from the viewpoint of a posteriori error estimation, it 
has been shown [10] that a type of equilibration condition has a significant impact on the 
efficiency of a posteriori error estimators for more general second order differential operators. 
In order to illustrate the idea let us consider the simple example: 

Find w such that 


— Aw(x) + w(x) = f{x) in fl subject to — — = g on d£l (48) 

The issue here is not the existence of solutions to the problem since the problem above is 
always well posed. However, when one is attempting to compute a posteriori error estimators 
for the problem, once again a local problem on each element K of the same form arises with 
data related to the residual / + Aw — w and the jumps in the normal fluxes across the 
interelement boundaries. In [10] it is shown that the resulting estimator can be rather 
pessimistic if the data for the local problem is not chosen with some care. Quantitatively, 
this statement means that the data gx for the local problem should be chosen such that 

J (/(*) - w{x)jdx + j>^g K {s)ds = 0 (49) 
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or at any rate such that this quantity is sufficiently small. If this were the only constraint on 
g K then there would be little difficulty. However, there is an additional constraint between 
g K and the flux g L on the shared side Y K l of the neighboring element which destroys the 
local nature of the condition. 

Once again, the problem is to split the fluxes on the interelement side between the two 
neighboring elements in such a way as to satisfy a type of generalized equilibration condition. 

Flux Recovery Scheme 

In this section we give a heuristic alternative interpretation of the flux splitting discussed 
above. Flux recovery and postprocessing schemes have been intensively studied [13] for the 
finite element method on uniform meshes for problems with smooth solutions. 

The strategy used in this work to obtain enhanced approximations to the flux was to 
average the fluxes along the interelement edges. Since the meshes were assumed to be 
essentially uniform and of uniform degree by symmetry, one should take the average of the 
two possible flux approximations along the edge. Indeed, one often finds that this type of 
scheme yields superconvergent approximations. 

For nonuniform meshes, however, a weighted average should be chosen instead of a simple 
average. Also, in order to effectively deal with non-smooth functions, the weights must be 
allowed to depend on the approximation itself, giving a nonlinear postprocessing. 

The question then remains as to how these weights should be chosen. Since the true 
solution satisfies (48), integrating over an element K gives 


[ ( rt \ 

JkV™ 


/ .. \\ .j 

j j u<k 


j dw 
fdi< dux 


J _ Cl 
AJLO — U 




Therefore, one might try to choose the approximation to the flux gx in such a way that 


(/(*) - w{xfjdx + j>^g K {s)ds = 0 (51) 

which is the same condition as was arrived at in the previous section. 

Notation and Preliminaries 

In order to simplify the presentation, we shall consider the two-dimensional case. The 
extension to three dimensions will become clear after a few remarks of clarification. 

Let Cl CM 2 be an open bounded domain with boundary T consisting of a finite number 
of smooth arcs meeting with internal angle 6 € (0,2tt). H m (Cl), m € Z + will be used to 
denote the standard Sobolev space equipped with the norm ||u|| m> n. We use the notation 
H°(Cl) = L 2 {Ci) in the case m = 0. 

Let V denote a partitioning of Cl into a collection oi N = N(V) subdomains K such that 


1. N(V) < oo 
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2. ft = (J I< 

— * 

3. If K / Z, than I< n L is empty 

4. I\ are Lipschitzian domains with piecewise smooth boundaries dK 
5- F KL = dK D dL is an entire edge of at least one of I< or L 

It is convenient to formally define the exterior of 0 to be the zeroth element. In this way, 
the complete set of element edges may be characterized by 

E = E{ V) = U ( 52 ) 

~K,L-.K>L> 0 

The set of interelement edges may be characterized by 

E I = E I {V)= U t kl ( 53 ) 

-+K,L-.K>L> 0 

The unit outward pointing normal vector on I< is denoted by n K . Let 


j +1, if I< > L, 
\ -1, if I\ < L 


and define 


(54) 


n(s) = <rKLni<{ s ) = ^lk t ^l( s )i s £ 

That is, n points outward from the domain with the largest index. 

We now focus our attention on a generalized model problem. Let n be the number of 
components of the problem and x = [# W]"- Let B: X x X -» * denote the bilinear form 


B(u,v) = ^ B :k (uj,v k ) 
},k=l 


(56) 


where u = («i, . . . , u«), v = (t>i» • • • > v n) an d 


B jk (uj, v k ) = {Vujt • A kj Vu 3 + v k B k] ■ Vuj + v k C kj Uj + Vv k ■ D kj u 3 ) dx 


+ 


L{ v ' b -ih +vkct ’ u ’} ds 


(57) 


22 


where A k > € R 2 * 2 , B kj € R 2x \ D k > € &'* 2 , and C kj , b k \ and c k J € M. We shall assume 
that the bilinear form is continuous and coercive on x x AS that is to say, there exist positive 
constants M and a such that 


B(u,v) < M ||u|| x ||r|| x V € x ( 58 ) 

and 


B(v,v) > a||v||* V r 6 x 
Let £: x — » R denote the linear form 

£( v ) = E 

it=i 

where 


(59) 


(60) 


L k (vt) = [ f k Vkdx + [ f k v k ds (61) 

J 0 JdQ 

where f k and f k eR. We shall assume that the linear form is continuous on x- 

Under the above assumptions there exists a unique solution u 6 x of the problem: 

Find u G x suc h that 


B{u,v) = C(v) V v G x (62) 

Let x C x denote a finite dimensional subspace consisting of continuous piecewise polynomial 
functions defined on the partition V. The polynomial degree is allowed to vary from element 
to element but the functions are constrained in such a way that interelement continuity is 

preserved. 

Our error estimation problem consists of estimating the error in an approximate solution 
w G X where u is obtained using some approximation technique. It will be useful to consider 
the structure of the space x in ore carefully. In particular, let A be any unconstrained or 
proper vertex A in the partition. Associated with each such node is a piecewise polynomial 
of degree one, which vanishes at every other regular node in the partition and (subject 
to suitable scaling) has the value 1 at the node A. For fc-irregular meshes [14], the support 
of ip A consists of a patch of elements containing or near to the node A. In the case of regular 
meshes, the support is merely the set of elements of which A is a vertex. 

Let F(V) denote the set of all unconstrained nodes in the partition. It is readily seen 
that the set { i{> a }a er forms a partition of unity subordinate to the covering {I<} on the 
domain Cl. That is to say 
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(63) 


E M*) = h Vxen 

AeT(T) 

Let us note that (62) implies a fortiori 

Y,B’ l (u,M) = L k (M < 64 > 

for fc = 1, • . . , n and any A G ^(V). 

For each element A G 'P and for j = 1 , . . • , n let 

Ct (*) = |£ (^'Vu; + D kl Uj) J (65) 

i.e., Q k K denotes the fc-th component of the flux on element K. With each interelement edge 
Tkl € E(V) we associate n functions ajf[(s), 5 G T/vtl corresponding to each component of 
the flux. We shall distinguish between the and ct^', and usually these will be unequal, 
but will be required to satisfy the condition 

«!&(») + = !. s ^ r Ki- <66) 

The normal component of the flux will usually be discontinuous across the interelement 
edges. The jump in the flux is denoted by 

[n • Q*]] (s) = n K ■ (Q k K - Q k L ) = n K ■ Q k K + n L ■ Q k L , s T K l (67) 

The a^ L are used to construct an average normal flux along the T KL from Q k K and Q k L as 
follows 

(n ■ Q k ) 1 _ a = {blkQk + q klQl ) * n 

Remark: The subscript is written as 1 - a to emphasize the fact that the average is formed 
using a LK times the flux from element K rather than the more natural notational choice of 

<*KL- 

Existence of Self-Equilibrating Flux Splittings 

The preceding remarks indicate some of the reasons for obtaining self-equilibrating flux 
splittings. However, the question of whether such splittings exist arises. In this section we 
prove that such splittings can be found. 

First of all, we formulate the appropriate equilibration condition for the model system 
(62). The strong form of (62) is to find («i, . . . , u n ) such that for k = 1 , . . . , n, 
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{-V • (A kj Vuj) + B kj • V Uj + C kj Uj - V • (D^Uj)} = f k in 0 (69) 

i=i 

and subject to the boundary conditions for k = 1, . . . , n 

£ |n • (A fcj Vuj + ^ + c kj Uj | = /* on dfl (70) 

Let I< eV and let <j> be sufficiently smooth on K, then multiplying (69) by (f>, integrating 
over I< and applying the boundary condition (70) on dl< R dfl, gives for k = 1, . . . ,n 

± B% ^A) = 4 (« + L m ^ K ■ £ + *>*“>) ds < 71 > 

J=1 J oh \au J_! 

and now choosing (j) = 1 gives 

£ Bjf 1) = 4(1) + L, »K • E ds < 72) 

j=i j=1 

This equation expresses the equilibration of the true solution to the model problem. The 
equilibration condition for approximation u and hence for the error estimation problem is 
correspondingly 


£b*(S j ,1) = £1(1) + l KV „( nK Qt \-« d3 (73) 

i— 1 

Although the above derivation is rather informal, the result can be obtained rigorously by 
following the corresponding steps for the scalar case in [10]. 

Let denote the lack of equilibration in the k - th component on element K, that is 


Ak l) = ltd) - £ Bi + / 8KVM <»K • *• ( 74 > 

The remainder of this section is concerned with proving that there exist choices of a for which 
the equilibration condition can be satisfied on every element K € V and all k = 1, . . . ,n. 

The replacement of unity in (74) with the partition given by (63) yields 


where 




' K 


= E A 


K,A 


A 


(75) 
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Therefore, a sufficient condition for equilibration is to ensure that vanishes for all K , A , 
and k. Rather than work with a$ L directly, it will be convenient to introduce the following 
quantities 

(k) ( k ) 1 i (fc) (k) ^ (77\ 

Ph'L = a l<L ~ 2 and PLK = Ct LK- g V 77 ) 

so that using (66), we obtain that n is anti-symmetric: 


{k) ( k ) 

Ml = -nlc 


Using the new variables, one finds that 

(% • Q k ) j_ o = ( n K ■ Q k )i_ + PLK [[ n • 


and that 


(78) 


(79) 




j= 1 


/ M s ) (n K • Q k ), ds + [ Ms)Pkl [[ n A’ * Q 

JdK\dQ ' '2 L>o Jt kl lL 


ds 


Denote 


and 


Plk,a = “ / t a( s ) [[ n A- * Q k ]\ ds 

J 'KL 

then the vanishing of is equivalent to 

£ = iZ 

L> 0 

Let = PlkPlk,a then we obtain 


(80) 


^ k]a — (Vu) — XI {Uji^A) + / ^aIs) ( n K Q )i ds ( 81 ) 

Jdl\\du 2 


(82) 


(83) 


V U {k) - b {k) 

/ , PLK — °K,A 


L> 0 


and using (78) to eliminate for L > K gives 


( 84 ) 
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( 85 ) 


E ( ~(k) ^(k) \ ilk) 

(PLK “ Vkl) = b K,A 

K>L> 0 

which may be written in the form of a sequence of linear systems of equations 

M A fiM = b (k) (86) 

for k = The underlying Matrix M A depends only on the topology of the patch 

of elements forming the support of xp A , and not on the index k. We shall return to this 
point later since it has a significant impact on the cost of the algorithm for computing the 
splittings. 

The structure of the matrix M A was studied in [10] where it was shown that the kernel 
of the adjoint matrix M m A was the one-dimensional subspace spanned by vectors of the form 
(1,...,1). Thus a necessary and sufficient condition for the existence of solutions to the 
system (86) is that for k = 1, . . . , n: 


£ C. = o 

K$P 

Now for the special case where u satisfies 

B(u, v) = C(v) V v € x 

then 


(87) 


(88) 


EC, = E 

KeP KeP l. j =1 ' * 

= IHM - £ Bit ft.M+L L. m Ms) ■ Q k )i & 

j= i K& ,Jdh\dn 2 

= Yi I Ms)(n K -Q k ) i ds 
fa, JdK\dO ' ' 5 

= H / '*I > a{. s ) (T kl ( n • Q k )i d s 

KePL€P JrKL 2 


(89) 


= 0 

where we have made use of (68) and (54). Consequently, for this case, the systems (86) 
always have a solution which is unique up to the addition of an arbitrary multiple of the 
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vector (1, . . . , 1). Consequently, we have shown that can be made to vanish and thus 

that there exist choices of o; for which the equilibration condition is satisfied. 

An Algorithm for Flux Splitting 

The principal results of the previous section are that the self-equilibrating flux splittings 
can be constructed directly using solutions of the systems (86). The matrix Ma does not 
depend on a particular component of the system, but only on the topology of the patch 
formed by the elements constituting the support of the linear basis function ip a- Therefore, 
for the purposes of this section we can dispense with the superscript k referring to the 
component of the system. 

Without loss of generality, we suppose the elements forming the patch to be numbered 
from 1 to Na- The matrix M A also depends on the interelement edges between the elements 
forming the patch. Let us arbitrarily label these edges from 1 to Ne- Then, Ma has one 
row corresponding to each of the elements 1, . . . , N A and one column for each of the edges 
T u ...,T N e . Examining the equations leading to the linear system (86), we see that the 

structure of Ma is 



0 

1 

-1 


if dSli fl dUk is empty for all k, 
if dfli H dflk = Tj and k < i 
if di 1, fl dflk = Tj and k > i 


The matrix Ma is singular, meaning that the solution p of 


(90) 


M A p = b A ( 91 ) 

is unique only up to the addition of a multiple of (1, . . . , 1). From a theoretical point of view 
it is irrelevant which solution of this system we pick. However, from a practical viewpoint, 
one should choose the smallest solution since this will inhibit the effect of rounding errors. 
Therefore, let us try to find p such that p T p is minimized subject to (91). Let y denote the 
Lagrangian 

y{p,\) = 7j£ T /l - {M A p - b A ) T A (92) 

then the conditions for a stationary point are 


p = MjA (93) 

and 

M A p = b A (94) 
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Together these conditions imply that 


M a M t a \ = b A (95) 

Therefore we may obtain p by solving (95) and then using (94) to compute p. At first sight 
this process seems to offer little over attempting to solve (91) directly. However, the special 
structure of the matrix T A defined by 


T a = M a M t a (96) 
means that there are advantages in pursuing this approach. 

Structure of T A 

First, notice that T A is symmetric and positive semi-definite. The kernel of T A is identical 
to the kernel of M A . Making use of (90) we readily conclude that the elements of T A are 
given by 


Ci if j = 1 


Mii = 


-1 

0 


if fl, and f Ik share a common edge, 
otherwise 


(97) 


where C, is the number of elements in the patch which share an edge with element fl,. It 
is apparent that T A can be rapidly constructed purely from the topology of the patch of 
elements forming the support of tp A , and for this reason we often refer to T A as the topology 
matrix for the node A. Some examples of topology matrices for various types of patches are 
shown in Figs. 3.2 and 3.3. 
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Solution of T A fi = b A 

The matrix T A is still singular and the same considerations mentioned earlier still apply. 
Moreover, since u may not satisfy (88) or due to rounding error and the familiar variational 
crimes [14], the right hand sides may not in fact satisfy the condition (87) precisely, 
meaning that the system actually seen by the computer is singular and has no solutions. 

In order to combat this type of difficulty we apply a two-step procedure: First, in order 
to ensure that the matrix equation has a solution, we modify b ^ according to 

C - C - (98) 

for I< = 1 , . . . , N A , where 

' A l> = 4-E b K,A ( 99 ) 

^A K=l 

This process means that the condition (87) is enforced while avoiding destroying too much 
of the information present in the components of b^. Second, in order to make the matrix 
equation easier to solve and to obtain a reasonably small solution from the set of solution 
vectors, we modify T A to become 

t a <-t a + d (100) 

where 


D = 


1 1 ••• 1 

1 1 ••• 1 


1 1 ••• 1 


( 101 ) 


This does not change the essential properties of the matrix owing to the condition (87) 
satisfied by the data, but merely picks a particular solution. However, it does mean that 
the matrix equation is easier to solve since it is now nonsingular. An alternative method for 
picking out a particular solution would be to set, for example (T ^4 )ix to some large value. This 
is similar to the “big spring” method used in finite element analysis. However, this would 
completely destroy the information present in the first component of b A \ The method chosen 
balances the information which is lost between all the components of b A \ 

To see that the matrix is nonsingular, let x (E R Na denote any vector. Let V = 
span(l, . . . , 1) and V 1 be the orthogonal complement with respect to the Euclidian inner 
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product on R Na . This * may be written as y + z where y 6 V is the orthogonal projection 
of * onto V and z = x - y G V x . Using these definitions it follows that 


x t (T a + D)x 
= x T (T A + D)(y + z) 

= x T (T A z + Dy) (102) 

= z T T A z + y T Dy 

> 0 

where we have used the fact that T A and D are symmetric and that V is the kernel of T A . 
Now since z t T a z and y T Dy vanish only if z and y are zero respectively, we see that the 
matrix T A + D is positive definite and hence nonsingular. Moreover, the condition number 
of the matrix is of order unity. Thus we compute the solution A of the system 

(T a + D)\ = b A (103) 

The flux splittings are obtained to be 


1 , PKL 
q kl = 77 + “ — 

2 Pkl 


(104) 


where fr is obtained from A using equation (93). Owing to the special structure of Ma, this 
Kn cim'TklifiArl fnrt.Vipr t.n crive. 

lUlillUiUi VWiJli ^ wli**i/a***VN* * Q" 


1 , Ak - A i 

®KL = 77 d 

2 Pkl 

The final splitting factor to be used in (68) is given by 


(105) 


^) = T.4UMs),ser KL doe) 

>1 

Of course, most of the terms in this summation vanish due to ip A having non-zero values 
on a small number of edges. For example, in the case of regular meshes, only two terms 
in the summation are non-zero, namely those corresponding to the two nodes forming the 
endpoints of the edge Ykl • In the case of irregular meshes, the situation is more complicated 
with at most three non- zero terms appearing in the sum. 

Model Problem Example 
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The following equations were solved on the square domain fi — , l) x (o, 2 ) 

dudu n 

- eA u + x— - y -x- + xyu + v = U 

ox oy 

dvdv n 

- eAv + x— - y-x- - u - xyv = 0 

ox oy 

dw dw 

~ eAu> + x Tx ~ "aF = 

subject to Dirichlet boundary conditions coinciding with the known analytic solution 


( 107 ) 


u 0 = e ( x2 - y2_1 ) /2e 

n 0 = xye^~^ ( 108 ) 

w 0 = xy 


where £ was chosen to be 0.01. 


Estimate Formulation 

This problem was solved on a sequence of /ip-adapted finite element meshes and the errors 
in the approximations were estimated using a version of the element residual technique with 
flux balancing. Since the formulation of the element residual technique is neither obvious 
nor unique due to the unsymmetric nature of this problem, we provide a brief derivation. 
First, the variational forms for the governing equations can be written 


f du du \ , 

eVu • V<£ i + lx— - y-^ + xyu + vj <f>i 


B 2 {u,(f> 2 ) = f 

J W 

r \ ( dw dw\ 

(«.*■) = /„ [ eV “ ' + ~ 

V 4 > = € X = l-H 1 (O)J 

where u = (u,u, w). By denoting 


( dv dv \ , 

eVv • V (f> 2 + lx— - y— - u - xyv J <f> 2 


dn = o 


d fl = 0 
dU = 0 


(109) 


the problem may be stated: 
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B(u,v) = ^2 Bi(u,v) 

i - 1 


( 110 ) 
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Find u€x such that 


B{u,<t>) = 0 v</»ex ( m ) 

The approximate solution comes from a subspace of \ an< i f° r this particular example, can 
be formulated as 

Find u € x C x suc h th a t 


B = 0 V (f) (E x (H 2 ) 

The error e = (e u ,e„,e w ) in the corresponding approximation solution u = ( u,v,w ) can be 
written 


v <j> e x 

For the purposes of error estimation, however, we define another bilinear form, 

~B(u, v ) — J [eV u ■ Vv + u • ujdfl 

Our error estimate, e, is the solution to 


/ — / \ / jl\ 

■D (e, <pj = 0{e, qj) = -u v'u, qs) 


U JL r 
v 'f' <= A 


We use B to define an energy-like norm 


( 113 ) 


( 114 ) 


n i ca 


|||.||| ! = B(v) ( 116 ) 

and under standard assumptions it can be shown that the energy-like norm of the estimate 
is equivalent to the error in the H l norm 

lllelll ~ Ml, (117) 

Thus, the global error estimate problem can be written 
Find e G X such that 

B{e,<t>) = n{u,<t>) v <f> a x ( 118 ) 

where due to the absence of body forces, the residual, 


n{u,(})) = -B{u,<j)) 


( 119 ) 
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Finally, the error estimate is 


, = ]||e||| = v/B(e,e) (120) 

In order to reduce this global problem to element-by-element problems, the constraint of m- 
terelement continuity of the error is relaxed and then reapplied using Lagrange multipliers. 
The Lagrange multiplier is roughly identifiable with the interelement flux which is approx- 
imated by the equilibrated fluxes of the approximate solution. Thus, the error estimation 
problem for an element k can be written 

Ti k {e,(j>) = n k {u,4>) + J^^(n k -Q k ) x _ a <t>ds ( 121 ) 

where the flux term was discussed in this section and the error estimate is given by 

T] k = Median* = \fBk (e,e) (122) 

Typically this local problem is solved over some prescribed subset of x and frequently over a 
subset of x/x (since 7 = 0 V <p G x)- We denote this subset as M and if x is the set 
of piecewise polynomials of order < n, then M is the set of piecewise polynomials of order 
m where n < m < n + 1 . 

Results 

This advection diffusion problem was solved on a sequence of /ip-adapted meshes. For 
each mesh the approximate solution, error estimates, and the true errors were calculated. 
Both the estimate and the true error were evaluated in the energy-like norm defined earlier. 

The beginning and final meshes are shown in Figs. 3.4 and 3.7 while gray scale plots 
of the error estimate and true error are shown in Figs. 3.5 and 3.8 and Figs. 3.6 and 3.9, 
respectively. These plots are indicative of the overall behavior of this estimate. Namely, the 
global estimated error is fairly close to the true error (especially for fine meshes) while the 
distributions of the error estimates are for the most part quite accurate. 
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Derivation of Error Estimation for Incompressible Flow 

The procedure followed in deriving an error estimate for incompressible flows is very 
similar to the derivation of the estimate presented for the model problem example with two 
added complications. First, since the governing equations are nonlinear, the formulation 
of the error is also nonlinear. Unfortunately, this nonlinearity destroys the equivalence 
between the error estimate and the error. The consequences of this lack of equivalence are 
not fully understood. In addition, the error in the momentum equation and the error in the 
incompressibility constraint must be balanced. The procedure 
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D.O.F= 25 


Figure 3.4: Initial mesh used in the model problem example. Mesh consists of 19 bilinear 
elements. 
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MIN=0.0034797 
MAX=0.204856 
ERR OR =0.4037959 
D.O.F= 25 


Figure 3.5: Estimated error distribution on the initial mesh for the model problem example. 
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MIN=0.001933 
MAX=0. 1444008 
ERROR=0.203888 
D.O.F= 25 


Figure 3.6: Actual error distribution on the initial mesh for the model problem example. 
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Figure 3.7: Mesh used in the model problem example after four h-p adaptive passes. Mesh 
consists of 25 elements ranging from bilinear through sixth polynomial order. 
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MIN* 0.524E-05 
MAX=0. 0020944 
ERROR=0.0033857 
D.O.F= 340 


Figure 3.8: Estimated error distribution on the adapted mesh for the model problem example. 
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by which this was accomplished for Stokes’ flow was described in the previous section 
and a similar technique is used for incompressible flow. 

First, the steady state incompressible Navier-Stokes equations are written in the weak 
form 


a(u,v) + c(u,u,v) + b(v, P) = 0 
b{u,q) = 0 V(v,g)€x x ^ 

where x = [f/ 1 ^)] 2 an d 

M= | 9 GL 2 (0): J u qdtt = oj 

and 


(123) 

(124) 


a(u^v) = —J + (Vu) 7 ^ : + (X7v) T ^jdQ, 

c(w,u,v) = j (w-\?)u-vdQ (125) 


b(v,P) = - f P(V * v)dSl 

J o 


rr. i.: ^ 

xaiaiifc x 


X aiiu ivi 


C M the corresponding approximate solution (u, P) fc X x M satisfies 


a(w,«) + c(u,u,v) + b(v,P^j =0 
b(u,q) = 0 V (ii,9)6xxM 

The formulation for the error, however, differs due to the nonlinear term 

a(e,t>) + c(e,e,i>) + c(e,u, v) 

+c(w,e,t>) + b(v,E) = 

—a(u,v) — c(u,u,v) — b (v,P) 
b{e,q) = 

where 


(126) 


(127) 


e = u — u 

E = P - P 
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MIN= 0.659E-05 
MAX=0.0017857 
ERROR=0.0028596 
D.O.F= 340 


Figure 3.9: Actual error distribution on the adapted mesh for the model problem example. 
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Our error estimate ^e, Epj is defined as the solution to 


where 


a(e,r) + (e, v) = 
+ 
+ 

( E,q ) = b(e,q ) 


a(e,v) + c(e,e,v) 
c(u,e,v ) + c(e,u,v) 
b(v,E) 

v ( V,q ) e x X M 


Correspondingly, we 


(tt, v) = / u ■ v dQ 

Jn 

define a weighted H 1 norm 


(128) 


IIHH = y/a(u,u) + («,«) 

At this point we would like to state that the weighted norm of the estimate is equivalent 
to the H 1 norm of the error but unfortunately this relationship is lost due to the nonlinearity. 

For a sufficiently accurate approximation, however, using u ~ u and e <C u and we can 
neglect the nonlinear term and recover the equivalence relationship. With this in mind, the 
error estimate is formulated 


a(e,v) + (e,v) 

(£,«) 


7l e ( U,P,V ) 
TZ e (u,q) 


where the resjduals are defined 


(129) 


1l e (u,P,v ) = -a{u,v)-c(u,u,v)-b(v,P) ^ 130 ^ 

1 lE(u,q) = ~b(u,q) 

As with the previous example, these global equations are reduced to element-by-element 
formulations and solved over some subspace of x/x- The resulting estimates are evaluated 
in terms of their corresponding norms 




k 

e 



n k 


Ve 



L 2(fi*) 


(131) 
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Using balancing arguments, these estimates are combined to yield the final error estimate 
for an element 



where e > 0 is a user— specified parameter weighting the two error estimates. For our partic- 
ular example, e was assigned the value unity. 

Results — Laminar Flow Over a Backstep 

This method of error estimation was applied to the solution of laminar flow over a back- 
step (Re = 200) on a mesh of 196 biquadratic elements. The stream function corresponding 
to this solution is shown in Fig. 3.10 while the estimated errors are shown in Fig. 3.11. For 
the purposes of comparison, a solution was computed on a mesh of 196 fifth order elements 
(having over 5000 degrees of freedom) and the difference between this highly accurate solu- 
tion and the solution on the original quadratic mesh was measured in terms of the weighted 
H 1 norm of the velocity difference and the L 2 norm of the pressure difference. The resulting 
plot is shown in Fig. 3.12. Note that the global quantities are within a factor of two of each 
other and the distributions of error estimates and solution differences are very similar. 
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Figure 3.10: Stream function corresponding to solution of laminar flow over a backstep at 
Re = 200 on a mesh of 196 biquadratic elements. 


47 



ERRORS 
MIN= 0.142E-03 
MAX=0. 02 07302 
GLOBAL=0.0822 
D.O.F= 849 


Figure 3.11: Error estimates for solution to laminar flow over a backstep at Re — 200 on a 
mesh of 196 biquadratic elements. 
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3.2.4 Error Estimation and Postprocessing for Multiple Problem Types 

In this section we identify the issues which arise when one attempts to construct a general 
purpose package which is capable of analyzing the solution obtained when an elliptic bound- 
ary value is approximated using a h-p finite element discretization. The aim is to provide 
techniques which are applicable to a very broad class of problems: 

• self-adjoint elliptic systems of partial differential equations (e.g. Lame-Navier equations 
of linear elasticity) 

• problems subject to a constraint (e.g. Stokes’ approximation to the incompressible 
Navier Stokes equations) 

• non-self-adjoint and convection dominated problems such as those arising in heat trans- 
fer and fluid flow 

• problems with strong convection and a constraint (e.g. Oseen’s approximation to the 
incompressible Navier Stokes equations) 

• non-linear problems with some or all of the above features (e.g. the incompressible 
Navier Stokes equations themselves). 

The main purpose of the techniques may be divided into two categories. 

1. To provide some numerical estimate of the discretization error or quality of the solution 
(i.e. to obtain a posteriori error estimates) 

2. To provide enhanced approximations over the original approximation to various derived 
quantities associated with the solution (e.g. to postprocess the solution). 

Moreover, the techniques should be able to cope with and perform effectively for general h-p 
finite element approximations. That is, the underlying mesh used in the discretization may 
contain irregular or constrained nodes and have elements on which the polynomial degree may 
vary from element to element. Additionally, the solution may consist of components which 
vary in degree of approximation and smoothness. For example, in the case of approximation 
of the Stokes equations one must allow for the pressure approximation to be discontinuous 
between the elements and also of lower order than the approximation to the velocities. 

Lastly, but by no means least, the technique must be reasonably economical. That is to say, 
the cost of the postprocessing should not exceed the cost of resolving the original problem 
using a refined discretization. 

The first question to be resolved is that of how the error should be assessed (i.e. in what 
norm it should be measured). Traditionally, the energy norm has always been chosen, mainly 
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Figure 3.12: Norms of differences between solutions to laminar flow over a backstep at 
Re = 200 from meshes of 196 biquadratic and 196 fifth order elements. 
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because of the difficulties of working with other norms. While it is convenient to work with 
the energy norm in the case of self-adjoint problems, for other classes of problem the choice 
is less clear cut since there is often no natural norm with which to work. Moreover, in the 
case of problems for which there is a constraint, one has to come to a decision on the correct 
way in which to deal with this feature. 

The question still remains of how the error should be estimated. Currently, the techniques 
used may be divided into two categories: those based on solving an error residual problem 
c.f. [10,9] and those based on a postprocessed solution c.f. [15,16]. The approaches which are 
based on using a postprocesses solution appear at first to be very attractive. Unfortunately, 
it has been found both in theory and in practice that the techniques perform poorly when 
either the mesh or the true solution is not sufficiently (very!) smooth. Moreover, as we shall 
see below, for the problem classes which we have identified above, the correct approach to 
postprocessing is not at all clear. Conversely, the error residual technique seems to perform 
reliably even in the cases when the mesh and/or solution are not smooth. However, it is not 
immediately clear how the approach can be extended to all of the classes we have discussed 
above. 

In summary, the following issues must be addressed when designing a general purpose error 
estimation package: 

• what norm can be used in the case of non-self-adjoint problems 

• how should the incompressibility constraint (e.g. in Stokes problem) be dealt with 

• what basic approach should be taken (i.e. error residual calculation versus postpro- 
cessing). 

In this section we shall give a rather complete answer to all of these problems associated 
with error estimation based on error residual methods. This is particularly pleasing from the 
point of view of developing a general purpose accuracy assessment package since it reduces 
the number of basic algorithms which must be present in the code. 

The issues associated with solution enhancement are less clear cut than in the case of error 
estimation. In part, this arises owing to the varied nature of the quantities in which the user 
is interested. All the user requires of an error estimator is that it provide an estimate of the 
global error and the local error distribution. 

The case of postprocessing is more complicated. For example, while performing an analysis 
of a linear elastic structure the user may require information on the strain or the stress at 
a point. Equally well, a derived quantity such as the stress intensity factor near a corner 
might be of interest. In solving a fluid flow problem, the user is interested in different 
types of quantities - in this case, primary quantities such as the velocity or pressure profiles. 
These issues have been addressed to some extent in the literature. For certain types of 
problems, such as linear elasticity, very effective techniques have been proposed based on 
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the use of extraction formulas [17,18,19]. Meanwhile, for the other classes of problem such 
techniques are as yet unavailable. It is noteworthy that in neither of the surveys [20,21] is 
there mention of recovery techniques suitable for the classes of problems described earlier. 
However, the preprint [22] does contain some ideas specifically tailored to fluid flow problems, 
but unfortunately does not readily extend to full h-p finite element discretizations. 

It can be anticipated that no single technique will provide the best application for every 
problem class. Thus, one needs to have flexibility to allow the user to input information 
according to the type of problem being solved and to the type of quantity to be recovered. 
The actual algorithm to be used would then be tailored based on both of these pieces of 
information. 

In the following sections we deal with each problem class in turn and outline how one can 
construct effective and reliable a posteriori error estimates and also discuss possible solution 
enhancement techniques. Numerical examples are given illustrating the behavior of the 
algorithms for various test problems. 

Linear Elastostatic Problems 

In this section, we discuss a posteriori error estimation for linear elastostatic problems. The 
purpose of this is twofold; first, to show how to obtain accurate a posteriori error estimates 
for h-p approximations to such problems, and; second, the theory developed in this section 
will be crucial in the subsequent sections. The material contained in this section is based on 
the work presented in [23]. 

Let us begin by considering a linearly-elastic body occupying an open bounded region fl C M 2 
with boundary dCl. The equations of equilibrium are written in the form 


-D'EDu 

u 

H{n)ED u 


/ 

in 

n 

0 

on 

r D 

9 

on 



(133) 


where D is the divergence operator, E is the elasticity matrix and H is the matrix formed 
using the components of the unit outward pointing normal vector n = (ni,n 2 ) as follows 


H(n) 


n-i 0 n 2 
0 n 2 


and g,f are the applied tractions and body force. 
It is convenient to introduce the bilinear form 


(134) 


B(u,v) = f (Du) l E(Dv)dx 
j n 

and the work done by the external forces is defined by the functional 


(135) 
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(136) 


L(v ) = J fvdx + g l vds 

With this notation the problem (133) can be rewritten as: 
find u € X such that for all v £ X 


B(u,v) = L(v) (137) 

The energy norm for this problem is defined by 

||u|| e=}/B{u,u) (138) 

The finite element approximation to u is obtained by performing a finite element discretiza- 
tion X of the space X and finding it 6 A such that for all v E X 


B(u,v) = L(v) (139) 

The error is defined by e = u - v. The goal is to construct an a posteriori error estimates 
of the error measured in the energy norm. Choosing v G X and using (137) gives 

B(e, v) = L(v) — B(u,v) (140) 

The quantity appearing in the right hand side of (140) contains the residual due to the 
approximation u fa u in the equations (133). 

Local Error Residual Problems 

Following the steps in the previous section, we reduce the single global problem stated in 
(133) into a sequence of local problems. As explained earlier, it is important to perform 
this decomposition in an appropriate manner if we are to obtain realistic error estimators. 
Specifically, it was shown how sensitive the effectiveness of the resulting error estimator 
will be to the way in which the interelement fluxes are decoupled. With this in mind, we 
introduce a quantity representing an approximation to the normal flux across an interelement 
boundary as follows 


(< 7 n (u)) « H(n)a (141) 

In the next section, we discuss various possible choices of interelement flux. With these 
considerations in mind, we now decouple (140) in error residual problems over each element 
K in the mesh: 

find <j>K £ Afft- such that for all w 6 Xk 

Bk(4> K ,v>) = l k( w ) ~ Bk{u,*>) + / w l {a n {u))ds (142) 
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Mesh 

True 

Error 

With Equilibration 

Without Equilibration 

Error 

Effectivity 

Error 

Effectivity 

1 

0.37113 

0.39692 

1.06949 

0.40592 

1.09373 

2 

0.24601 

0.26576 

1.08028 

0.26792 

1.08909 

3 

0.20319 

0.22043 

1.08486 

0.21823 

1.07405 

4 

0.12711 

0.16047 

1.26242 

0.14935 

1.17498 

5 

0.96581 (-1) 

0.12593 

1.30384 

0.10372 

1.07397 

6 

0.56591 (-1) 

0.88939(-l) 

1.57161 

0.63038(-l) 

1.11392 


Table lc Performance of Error Estimators for short cantilever example 


The notation B K (-, •) and L K { ) refers to the bilinear and linear forms defined in (135) and 
(136) evaluated over the subdomain I\ rather than over the whole domain ft, while the 
space Xk refers to the subspace consisting of the restrictions of functions in A to the single 
subdomain K. 

The approach is now to approximate the solution of (142) using a purely local finite element 
discretization as described in [24]. Having computed <j> K , the error on the element K is 
estimated as tx where 


e 2 K = B K (<l>K^K) ( 143 ) 

Approximation of Interelement Flux 

In order to determine an approximation of the interelement flux, we shall consider two 
approaches. It is worthwhile to bear in mind that along an interelement edge the finite 
element approximation of the normal flux is usually discontinuous owing to the discontinuity 
in the strains. Therefore, along any given edge, there are two possible fluxes corresponding 
to evaluating the flux using either of the two elements. One approach is to argue that a 
suitable approximation is to simply average the two contributions along an edge between 
the elements forming the edge. An alternative approach is to try to perform a weighted 
averaging of the fluxes and hopefully achieve an improvement over the simple averaging. 
Further details concerning the second approach will be found in [24,25]. In particular, for 
the latter scheme is is possible to show [10], that the error estimator which one obtains 
provides a guaranteed upper bound on the true discretization error. 

Numerical Examples 

In order to test the behavior of the error estimators described above, we apply each one to 
various test problems. 

Example: Short Cantilever 

The geometry of the problem is shown in Figure 3.13. 
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Mesh 

True 

Error 

With Equilibration 

Without Equilibration 

Error 

Effectivity 

Error 

Effectivity 

1 

0.53505 

0.60633 

1.13323 

0.59951 

1.12010 

2 

0.31535 

0.41129 

1.30420 

0.34988 

1.10949 

3 

0.28700 

0.38627 

1.34588 

0.32205 

1.12213 

4 

0.77801(-1) 

0.11322 

1.45510 

0.12039 

1.54723 


Table 2: Performance of Error Estimators for cracked panel example 


Figures 3.14 - 3.19 show the plots of the meshes used for computation and the resulting 
distributions of error estimates, while Table 1 lists the true and estimated (global) errors. 
Observe that as predicted theoretically, the estimates do provide an upper bound on the 
error. 

Example: Cracked Panel 

The geometry of the problem is shown in Figure 3.20 

This time we assume plane strain conditions and Poisson’s ratio of 0.3. The results of the 
adaptive analysis are shown in Table 2, while the sequence of meshes and corresponding 
error distributions are shown in Figures 3.21 - 3.24. Due to the symmetry of the problem 
it is only necessary to solve for half of the panel, which is loaded so that there is a stress 
singularity of the form r - 5 near the crack tip. As with the previous example, the results 
confirm that the error estimate provides an upper bound on the actual discretization error. 

Stokes’ and Oseen’s Equations 

In this section we consider the following model incompressible problems: 

• Stokes equations 


— i/A u + Vp = / 
V • u = 0 


(144) 


• Oseen’s equations 


— i/Au + (6 • V)u + Vp = / 

Vu = 0 (145) 

subject to appropriate boundary conditions. It will be convenient to write each of 
these problems in a generic weak form as follows: find u. p such that for all v and q 
there holds 
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L{v) 

0 


(146) 


B(u,v) + b{v,p) = 
b{u,q) = 

Remark: Strictly speaking the above formulation leaves much to be desired in terms 
of mathematical rigor, but our purposes and for ease of explanation it will suffice. 

These problems are natural testing grounds on the way to producing an algorithm 
suitable for the incompressible Navier Stokes equations. The main features introduced 
when considering these problems are 

- The question of what norm is appropriate to estimate the error for problems with 
an incompressibility constraint. 

- The finite element approximation to the pressure p may be discontinuous and of 
a different polynomial degree to the approximation to the velocity u. 

- The issue of how one should deal with the incompressibility constraint. Previously 
we have seen that it is appropriate to solve an error residual problem for each 
component of the solution, but it is not obvious that this is the correct direction 
in the case of a problem containing a constraint. 

- The introduction of a convective term in Oseen’s equations renders the problem 
severely non-self adjoint. It is unclear not only in which norm one should estimate 
the error, but also whether the convective effect will cause the local distribution 
of the error estimator to be faithful to the actual error distribution. 

Symmetrized Norms 

In this section we consider the question of the choice of norm in which to measure the error. 

Letting ii and p denote the h-p finite element approximation, we define the error in the 

velocity and the pressure by 


e = u — it 

E = p-p (147) 

One obvious choice of norm is given by 

||(e,£)||. = (||<*. +||£||1 ! )> (148) 

Here H 1 and L 2 denote the standard Sobolev spaces. In essence, the above norm weights the 
error in the gradients of the velocity against the error in the pressure. While this weighting 
is appropriate mathematically, it suffers the drawback that it is not readily accessible. In 
short, while it is desirable to use the star norm (148) it is not computable. 
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An alternative approach is required and one which suggests itself is to extend the concept of 
a symmetrizer which was successfully used in [26]. The idea which we present below extends 
the approach in [24,27] to problems involving a constraint. 

The first step is to introduce a new pair of bilinear forms 

A{u,v) = I t/(uij + Uj ti )vijdx 

(P,9) = / P<l dx ( 149 ) 

and then to define the pair <f>, rp to be solutions of the problems 


A((j>.v) = B(e,v) + b(v, E) 

(xp,q) = b(e,q ) (l 50 ) 

and then to consider the quantity 

||(*,0)|| = (^(<M) + (6V>))* ( 151 ) 

The quantity defined in (151) is referred to as the error measured in the symmetrized norm. 
It is easy to suspect that this quantity is related to the true error and this in indeed the 
case. The following result shows the precise relationship. 

Proposition 1 Suppose that the bilinear forms appearing in ( 146 ) are continuous and 
satisfy an inf-sup condition. Then there are constants C\ > 0 and C 2 > 0 for which 

Ci||(60)ll<ll(e^)ll.^^ll(60)ll ( 152 ) 

This result shows that while we are unable to make use of the start norm (148) directly the 
symmetrized error norm (151) is just as good. The next step is to show how we can estimate 
the error in the symmetrized norm. 

Error Estimation in Symmetrized Norms 

Examining the (150) and (151) it is easily seen that the error in the symmetrized norm can 
be rewritten as 

11(60)11* = A(<{>, 4>) + / (div ufdx (153) 

J 0 

This shows that if we want to estimate the error in the symmetrized norm, then it is only 
necessary to estimate the term A((j>,<t>). Returning to (150) and using (146) and (147) we 
obtain 
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Mesh 

Error Estimate 

True 

Error 

h- 1 

P 

Without Balancing 

With Balancing 

4 

2 

1.901(-3) 

1.778(-3) 

1.722(-3) 

16 

2 

4.64(-4) 

4.49(-4) 

4.15(-3) 

4 

3 

5.232(-5) 

5.229(-5) 

5.226(-5) 

16 

3 

6.528(-6) 

6.528(-6) 

6.524(-6) 


Table 3: Performance of Error Estimators for smooth Stokes example 


A(4>,v) = B(e,v) + b(v,E) 

= L(v) - B(u,v) - b(v,p) (154) 

The right hand of (154) is actually computable and represents the residual in the momentum 
equations from the original problem. A moment’s reflection on this equation reveals that 
the problem determining (/> is in fact a linear elastic problem owing to the definition of the 
bilinear form A(-,-). Having reached this stage, one can now compare problem (154) with 
the analogous problem (142) from the previous section. That is to say, one proceeds exactly 
as before to define local error residual problems for 4>. These can then be solved in the same 
manner yielding the first term in (153) and consequently an estimate of the error. 

Reliability of the Error Estimator 

The foregoing manipulations show how an a posteriori error estimate can be obtained by 
solving local linear elastic problems similar to those in the previous section. The question 
which naturally arises is that of how accurate are the error estimates. The following result 
quantifies the behavior 

Proposition 2 Let <j> denote the solution of the local error residual problems derived from 
(3.22) and define the local error estimate ex on element K by 

4 = + [ (div ufdx (155) 

Jl\ 

Then 

II(*.£)II.< Z t{E4-} ! < 156 ) 

The result shows that the estimator should bound the true error. 

Numerical Examples 
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In this section, we present numerical examples illustrating the effectiveness of the approach 
outlined above. 

Example: Stokes problem with smooth solution 

As a first example we solve Stokes problem when the true solution is known explicitly. While 
this is hardly a realistic problem, it does allow us to compute the true error exactly and thus 
be able to assess the performance of the error estimator. The true solution was chosen to be 


u 

v 

P 


cosx coshy 

-x 2 — sinx sinhy 

2 

vy 


(157) 


and the Stokes equations were solved on the unit square with both inflow and traction 
boundary conditions. The results obtained are shown in Table 3. The two types of estimators 
shown in the table correspond to two strategies for obtaining local problems from the global 
statement (154) (see [25] for full details). The results obtained indicate that the error 
estimators perform reliably. 

Example: Stokes’ Approximation to Backstep Problem 

The second example is the familiar backstep channel problem modeled using the Stokes’ 
approximation (i/ = 0.01). The geometry and boundary data are shown in Figure 3.25. 

Of course we cannot compare the estimated error with the true error in this case because the 
true solution is unknown. This time we have used the error estimator to solve the problem 
adaptively giving the sequence of meshes and local error distributions shown in Figures 3.26 
through 3.34. Physical intuition indicates the satisfactory behavior of the error estimator. 

Example: Oseen’s Approximation to Driven Cavity Problem 

In order to study the effects of a convection dominated flow, the Oseen approximation 
[v = 0.01) to the driven cavity problem is analyzed. Here, a standard driven cavity problem 
has been taken and a rotating convective field b = (y, -x) has been superimposed (see Figure 
3.35). 

Once again we cannot compare the estimated error with the true error in this case because 
the true solution is unknown. The error estimator has been used to solve the problem 
adaptively giving the sequence of meshes and local error distributions shown in Figures 3.36 
through 3.43. The performance of the estimator once again appears satisfactory in spite of 
the strong convective effects present in the problem. 

Summary and Conclusions 

We have outlined a procedure whereby reliable a posteriori error estimates can be constructed 
for h-p finite element approximations of broad classes of problems. The basic approach is 
the same for all types of problems and therefore is suitable as a basis for a general quality 
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assessment package incorporating uniformity of approach with reliability and robustness. 

The situation regarding solution enhancement is much less obvious. This is reflected in 
the literature where work on postprocessing of solutions of many classes of problems is 
either scarce or non-existent. For some classes of problems (such as linear elastostatics) 
very effective procedures are available. However, the approaches are extremely problem 
dependent. 

Faced with the problem of developing a general solution enhancement package, it is necessary 
to choose whether to require that the user input the type of problem being solved and to 
implement those procedures which are know, or, alternatively to attempt to construct a 
general approach. The advantages of the latter are obvious, but these must be carefully 
weighed against the disadvantages. In particular, it is doubtful whether any single approach 
can perform as well as the tailor made algorithms and, indeed, may even perform extremely 
poorly. If this latter course is to be adopted then the most promising area seems to be based 
on the element residual method. 
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Figure 3.14: First mesh and plot of error estimates of resulting solution for the short can- 
tilever. 
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Figure 3.15: Second mesh and plot of error estimates of resulting solution for the short 
cantilever. 
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Figure 3.16: Third mesh and plot of error estimates of resulting solution for the short 
cantilever. 
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D.O.r= 113 
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Figure 3.17: 
cantilever. 


Fourth mesh and plot of error estimates of resulting solution for the short 


65 





MESH 



D.O.F- 169 


Error Estimator (With Equilibration) 


i 



0 0.006 0.0U 0.024 0.03 


MEN=0.001903* 
MAX =0.02 9098 
D.O.F= 169 


Figure 3.18: Fifth mesh and plot of error estimates of resulting solution for the short can- 
tilever. 
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MESH 



D.O 224 



MIN =0.0012551 
MAX =0.017050 
D.O.F= 224 


Figure 3.19: Sixth mesh and plot of error estimates of the resulting solution for the short 
cantilever. 
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Figure 3.20: Geometry for cracked panel problem. 
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Figure 3.21: First mesh and plot of error estimates of resulting solution for the cracked panel. 
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MESH 





MIN=0.0055812 
MAX=0. 172894 
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Figure 3.22: Second mesh and plot of error estimates of the resulting solution for the cracked 
panel. 
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Figure 3.23: Third mesh and plot of error estimates of the resulting solution for the cracked 
panel. 
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Figure 3.24: Fourth mesh and plot of error estimates of resulting solution for the cracked 
panel. 
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Figure 3.25: Geometry for the backstep problem. 
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PROJECT: backstep Error Estimator (Without Equilibration) 
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Figure 3.27: Plot of error estimates of solution on first mesh without using flux equilibration. 
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PROJECT: backslen Error Estimator (With Equilibration) 
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Figure 3.28: Plot of error estimates of solution on first mesh with using flux equilibration. 
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PROJECT: backstop - ME.SI I - 
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Figure 3.29: Second mesh for backstep problem. 
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PROJECT: backstep Error Estimator (Wilhoul Equilibration) 





Figure 3.30: Plot of error estimates of solution of second mesh without using flux equilibra- 
tion. 
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PROJECT: backstep Error Estimator (Willi Equilibration) 
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Figure 3.31: Plot of error estimates of solution on second mesh with using flux equilibration. 
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Figure 3.32: Third mesh for backstep problem. 
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PROJECT: backstep Error Eslimalor (Wilhout Equilibration) 
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Figure 3.33: Plot of error estimates of solution on third mesh without using flux equilibration. 
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PROJECT: backslep Error Estimator (With Equilibration) 
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Figure 3.34: Plot of error estimates of solution on third mesh with using flux equilibration. 
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Figure 3.35: Geometry for driven cavity problem. 
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PROJECT: cav - MESII - 



Figure 3.36: First mesh for driven cavity problem. 
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Figure 3.37: Plot of error estimates of solution on first mesh without using flux equilibration. 
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Figure 3.38: Second mesh for driven cavity problem. 
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Figure 3.39: Plot of error estimates of solution of second mesh without using flux equilibra- 
tion. 
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Figure 3.40: Plot of error estimates of solution on second mesh with using flux equilibration. 
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Figure 3.41: Third mesh for driven cavity problem. 
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Figure 3.42: Plot of error estimates of solution on third mesh without using flux equilibration. 
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3.2.5 Residual Error Estimation: A Summary 

The residual technique for estimating errors in approximate solutions was discussed in the 
previous sections. This method seeks to estimate the error e = u — u k in an approximate 
solution Uh, where the analytic solution u = (ui,u 2 ,... u n ) satisfies an elliptic system of 
partial differential equations, which in their general form can be written 



+t»« £;+<*-/ 


subject to the generalized boundary conditions 


(158) 


+j+& + ®-£ +c —'* <i59) 

where the tensors A ke , B t , C, B n , B u B s , C b , f and f b are known functions of x = (x, y, z) 
and possibly m, and depend upon the particular problem being solved. 

After multiplying by a vector of test functions v, integrating over the domain (or bound- 
ary), and integrating some of the terms by parts, the combined governing equations and 
boundary conditions can be written in their variational form as: 


B(u, v) = F(v) VveH'itl) = H^Sl) x ... x H\U) 

s N r- 

n times 


( 160 ) 


where: 


B(u,v) = J n 

L 


dv T du T _ du ^ dv T j ^ 

t^i dx k dx ( £i dxt k= i dx k 


dx (161) 


+ 

+ 


an 

dv T 


t « du 7 ’ _ du T r% - — 

vB -^ + v B 'm +V B ‘~ 

dv T 


du 

7te 


dv T _ _ 

. D n u + -5 — Dtu -f -z—D $ u 4- C b u 
dn dt os 


ds 


and 



(162) 
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The corresponding approximate problem used to estimate the error in a given cell or element 
is formulated: 


Find 4> h tMh (fit) such that 

— f 3 

B k (<j> h , w)~ Y wTd ( n i <t> ds 

Jd(i k \an 

= F k {w) - B k ( u h , w) — / 

J d 


(163) 


do k \dn 


l t = 1 


3 

r 

e=\ 


Y w T g e n t - Y w T D(Ti(Uh 


ds 


du h 


+ Y I <*ew T Aij^—iiids VweMhittk) 

^ JdU t n (dQ k \9Q) OXj 


where Bk and Fk are equations 161 and 162 integrated over flk (instead of f!) and dflk D d£l 
(instead of dfl). Bk is the same as Bk if Bk is linear, otherwise, Bk is an appropriate 
linearization of Bk ■ The last term represents a weighted average of the element’s and its 
neighbors’ boundary fluxes. Several possible choices exist for the weighting coefficients, cxf. 


• For no average: ot( = (a*, a^, . . . a*); ay = 6(k 

• For a straight average: ay = (ay, ay, . . .ay); ay = \ 

• For a balanced average: ay = (ayi, ay 2 , . . . ay n ); each ay, is a piecewise linear function 
of s obtained by a flux-balancing procedure (see the previous sections.) 


The approximate problem represented by equation 163 is solved over a space of poly- 
nomial “bubble” functions, which is constructed by first taking the approximation space 
containing the solution UhtHh, p and adding polynomial shape functions until some new com- 
plete spectral order is spanned. These added shape functions form the basis for functions in 
the space Mh • Thus, if the approximate solution consists of (piecewise) polynomials of order 
p with a corresponding space Hi 1<P and if the enriched space contains polynomials of order 
p+ 1 (denoted Hh, P + i), then the three spaces are related according to Hh tP + 1 = Hh, P @Mh- 
For example, if the approximate solution on a brick element is trilinear (8 linear shape func- 
tions) and the enriched space is triquadratic (8 linear, 1 fully quadratic and 18 partially 
quadratic/partially linear shape functions), then the bubble function space Mh is spanned 
by the 19 additional quadratic shape functions. 

Once the error indicator functions <f> have been calculated, the error estimate is computed 
as a norm of <f>. Namely, 


e„ = |||<^||| = 


(164) 
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where the precise form of £?(•,•) depends primarily upon the type of norm desired. For 
energy type norms, B (•, •) may be very similar to other forms, B (•, •) and B (•, •). In fact, 
for symmetric linear problems such as heat transfer or linear elasticity problems, all three 
forms are identical. 

In summary, the process of computing residual error estimates for a given mesh consists 
of four steps for each element. 

1. Compute the contributions of the domain integrals of equations 161 and 162 to the left 
and right hand sides of the linear algebra problem (henceforth denoted as “stiffness 
matrix” and “load vector”). 

2. Compute the contributions of the various boundary integrals of equations 161, 162, 
and 163 to the stiffness matrix and load vector. 

3. Solve the linear algebra problem to obtain the error indicators. 

4. Evaluate the norms of the indicators to get the final error estimates. 

5. Evaluate the mesh-wise norm of the solution to obtain a scaling factor. 

6. Normalize each error estimate by the scaling factor. 

3.3 Solution Enhancement 

This section describes the work done in the area of solution enhancement techniques. It 
includes summaries of the results of a literature survey and research performed on this topic. 

3.3.1 An Overview of Postprocessing Techniques 

One paradoxical feature of the finite element method is that it often delivers a rather poor 
approximation to the quantities of the most interest. For example, when solving an elasticity 
problem, we are primarily concerned with the stresses. However, the finite element method 
provides an approximation of the displacement. Obviously, by a process of differentiation 
one can use the approximation to the displacement to obtain an approximation to the stress. 
Unfortunately, this approach often results in a very poor approximation to the stress. In 
particular, the stresses will be discontinuous across the interelement boundaries. 

This characteristic of the finite element method has been extensively studied and conse- 
quently various techniques are often used to enhance stress approximations. Such techniques 
are usually referred to as postprocessing techniques. 

In this section we shall be concerned with attempting to find suitable postprocessing 
techniques for use with various classes of approximations including finite difference, finite 
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volume, and finite element approximations. As we shall see, many of the existing techniques 
are applied to very specific types of finite element approximation schemes and in some cases 
are strongly problem-dependent. It is not only that the schemes have not been analyzed 
mathematically— they just fail completely when applied to other situations. 

The remainder of this section is organized as follows. First, we give a survey of the 
most popular and successful postprocessing techniques in use today. Second, we propose 
several techniques which hold the greatest promise in the context of general approximations. 
Throughout this and the following sections we shall illustrate the basic concepts by using an 
example of finite element discretization of a simple model problem. 

Find u such that 


—A u + u = / 

in 

ft 

u = 0 

on 

To 

du 

dii ~ 9 

on 

r N 


(165) 


where ft is a domain in R 2 with boundary dft = T n U T D and where T n n T D is empty. We 
shall let u denote a standard finite element approximation to u. The actual details of the 
finite element approximation scheme will be made clear in each case, and will be assumed 
to be general h-p finite element approximation otherwise. 


Superconvergence Phenomenon 

Perhaps the most natural types of postprocessing schemes have their roots in the phe- 
nomenon known as superconvergence. Indeed, this phenomenon is discussed even in the early 
work of Strang and Fix [28]. Probably the work of Barlow [29] was one of the first to make 
use of the idea that at certain “exceptional” or so-called “stress” points, the accuracy of the 
direct approximation Vtt to the gradient is significantly higher than the accuracy at other 
points. In order to give a concrete illustration we consider the case of our model problem in 
one dimension. 


-u" + u = / on [0, 1] ^ 166 ^ 

u(0) = u(l) = 0 

We consider a uniform mesh of size h > 0 and uniform polynomial degree p on each element, 
see Fig. 3.44. 

Let u denote the finite element approximation. In order to postprocess the gradient u' 
on element f l k we evaluate u' at the mapped Gauss-Legendre points. That is to say, we let 
£i>& 2 , • •• ,£p denote the zeros of the p-th degree Legendre polynomial. 

These lie on the interval [-1,1] which we map to element ft fc using the transformation 
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(167) 


M-U] -» «*, m) = - 0 + jMi +0 

where ft* = [a, 6]. We denote = F fc (£,•)• Then the average error in the approximation to 
the true gradient at these mapped Gauss points decreases at a rate of an order of h greater 
than the approximation elsewhere. Expressed precisely, the result is 

l 

{ E x> K - «' \ x f ] ] ri - ChP+1 (168) 

U=u=i J 

whereas at other points the rate of convergence is 0{h v ). In short, u' is an order of h 
more accurate at the mapped Gauss points. 

Of course having identified the points at which superconvergence occurs, we can obtain 
a super accurate solution at any other point in the domain by interpolating through these 
special points. Figure 3.45 shows the scheme in the case p = 1. 

If we let G( u') denote the interpolated gradient, then it has been shown in [15] under 
certain assumptions: 

IK-G(«')ll wtl )<CMW'-"llMn) ( 169 > 

and that 

\\n"-G'(u f )\\ Lm <Ch p (170) 

This one-dimensional example illustrates the key ideas used by this technique. Unfortu- 
nately, their extension into higher dimensions is not straightforward and will be discussed 
for several two-dimensional element types in the next sections. 

Linear Triangles 

Consider the case of an approximation using polynomial degree 1 shape functions on 
triangular elements. It is well known [30,31] that there are no points at which the supercon- 
vergence is observed. Indeed, it was thought that it was impossible to obtain superconver- 
gence, until the work of Levine [30] revealed that a more subtle approach based on averaging 
gradients between neighboring elements can be employed. In order to explain Levine s idea, 
we imagine a pair of triangular elements in a large mesh, see Fig. 3.46. 

Levine showed that a superconvergent gradient can be obtained at the midpoint of the 
interelement boundary as follows: define 

G (Vu) (**) = i{Vu| ni (z*) + Vu| fi2 (a;*)} (171) 
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where Vu|n, denotes the finite element approximation to the gradient on element fii. It 
turns out that G?(Vu)(**) is a superconvergent approximation. The assumptions needed to 
show this result mathematically are that u is smooth (i.e., u € H 3 (Q)) and that the pair of 
triangles form a parallelogram. 

Translated into a statement regarding the whole mesh, the second assumption leads to 
the famous “six triangle property.” Namely, that each interior node in the mesh must be 
surrounded by precisely six triangles. Examples of meshes which satisfy and fail this property 
are shown in Fig. 3.47. 

One is severely tempted to dismiss the “six triangle” property as an artificial mathemat- 
ical requirement which can be safely ignored in practice. Unfortunately, it is necessary for 
superconvergence. Without this property one does not obtain superconvergence. 

Quadratic Basis Functions on Triangles 

The situation here is very similar to that for linear triangles. Namely, one is forced 
to average between pairs of triangles to obtain superconvergence. Correspondingly, the six 
triangle property must also hold. Figure 3.48 shows the stress points for quadratic elements. 

It has been shown [32] that the two Gauss points along the interelement edges can be 
used to obtain superconvergence. 

G (Vu) (a:*;) = i {Vu |n, (**,•) + Vu| n2 (**,•)} * = 1 > 2 ( 172 ) 

For the result to hold, one must assume that u is very smooth (i.e., u € H 4 (£l)). Once again, 
without the smoothness assumption and the six triangle property, the superconvergence is 

J U. J J 

UC8 bi uj'Cu ui acvci cij uc^taucu. 

High-Order Basis Functions on Linear Triangles 

For the case of cubic, or higher order basis functions on triangles, one is tempted to 
conjecture that the Gauss points on the interelement edge are the appropriate points where 
superconvergence occurs. Unfortunately this does not appear to be the case. To date, no one 
has been able to demonstrate superconvergence for cubic or higher order triangular elements. 
The general opinion in the “superconvergence fraternity” is that it does not occur and that 
it is impossible. 

Quadrilateral Elements 

While this superconvergence technique seems to be of limited usefulness for triangular 
elements, the same is not true for quadrilateral elements. The fundamental results for these 
types of elements were shown first by Lesaint and Zlamal [33]. Unlike triangles, superconver- 
gence occurs for all orders of quadrilateral elements. However, as with triangular elements, 
the solution u must be very regular (u € tf p+2 (fi)) and also the mesh must be very regular. 
Essentially the elements must be parallelograms and the overall partitioning must be topo- 
logically equivalent to a mesh of the form shown in Fig. 3.49. That is, there must be four 
quadrilaterals around each node implying a structured mesh. 
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Given a single element Q, k from the mesh and assuming that it is a piecewise polynomial 
of degree p , the points at which superconvergence occurs are at the images of the pairs of 
Gauss- Legendre points. That is, if F k 'Q —* flfc where Q = [—1,1] x [—LI] iS fhe standard 
master element, then the superconvergent stress points are at Fk{£i,tj) 1 — Li - P w ^ ere 
{£} are the zeros of the p-th degree Legendre polynomial. 

Summary of Superconvergence Phenomenon 

The primary advantages of using superconvergence as a basis for postprocessing are 
its ease of implementation and that it is extremely inexpensive. Also, just because the 
technique is inexpensive does not imply that it offers a limited improvement to the solution. 
When the assumptions are satisfied, the results are sometimes startling. However, the mam 
disadvantages are the regularity requirements on the mesh and the solution. 

Projection Methods 

While the superconvergence or averaging approach evolved largely from the mathematical 
community, an alternative technique has evolved from the engineering community. The basic 
idea behind the method can be explained as follows: 

We require a continuous approximation to the gradient, but Vu is discontinuous across 
interelement boundaries. Let • • ■ ,<Pn denote the shape functions used to obtain it. 

Then we write our approximation G(u) to the gradient in the form 

G(u)(x) = ( 173 ) 

where a u <* 2 , • • • ,<*jv and . . . ,Pn are constants to be determined. Since G will auto- 

matically be a continuous approximation to the gradient the question is how to determine 
these multipliers. A projection technique is used to determine 

J ( ttl , . . . , a N ] ft, . . . , M = jf {G{u){x) - Vti(x)} 2 dx (174) 

Writing down the conditions for a minimum point of J we obtain the equations 

Mot — g x (175) 

and 

Mfi=g y ( 176 ) 

where M is the N x N matrix with elements 

Mij = / tpi{x)<pj(x)dx (177) 

Jo 

and g x ,g y are the TV- vectors with elements 
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( 178 ) 


[^x]l — 


\9y]i = 

S a ^ x) % ix 


The matrix M is guaranteed to be non-singular and hence a nd /3 are uniquely determined 
from these equations. This approach is identical to the approach known as the “Conjugate 
Flux Method” of [34]. In the above notation and context, the method was described by 
Hinton and Campbell [35]. Surprisingly enough, there is little or no mathematical analysis 
to support this method [36]. However, numerical evidence shows that the approach can be 
extremely effective. A related method was presented and analyzed exhaustively by Bramble 
and Schatz [37], although they made very stringent requirements on the mesh without which 
it is not clear how their approach can even be implemented. 

More recently, this approach has gained in popularity owing to its application to error 
estimation using the Zienkiewicz-Zhu method [16,38]. 

Analysis presented in [36,37] along with extensive numerical experience suggests that the 
method does indeed yield superconvergence but not globally. Namely, near the boundary 
of the domain, the method does not only fail to enhance the accuracy, but even degrades 
it. Unfortunately, it is usually near the boundary where we are the most interested in the 
gradients or stresses, since this is where singularities occur in applications such as elasticity. 

This phenomena is illustrated in the context of error estimation in Figures 3.50 - 3.54 for 
a three-dimensional linear elasticity problem. This problem corresponds to the pure bending 
of a prismatic beam. Figure 3.50 shows a sketch of the beam and its loading configuration 
along with a picture of the deformed shape of the body. Since the analytic solution to this 
problem is a second-order polynomial, only meshes of linear elements were used. Figures 3.51 
- 3.54 show plots of the error estimates and analytic errors measured in the energy norm for 
a sequence of meshes. The energy corresponding to the analytic solution is 223.6068. Based 
upon these figures, we see that this technique provides reasonable estimates of the global 
error but poor indications of how that error is distributed. This is especially true in elements 
near the ends of the beam. 

A further disadvantage of the method is the vast expense entailed in solving the matrix 
equations for a and /3. In fact, for two-dimensional problems the cost of postprocessing is 
roughly twice the cost of obtaining vl itself. 

Cantin-Loubignac Iteration 

Let us consider the linear elasticity problem: 

—S l DSu = / in ft ' 

u = 0 on Tc > (179) 

Her = g on Tyv 
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where u is the displacement, <r = DSu is the stress, D is the elasticity matrix, S is the 
operator 


dx 0 
0 dy 
L dy dx J 

and H is formed from the unit normal vector n = (n x , n y ) as follows: 


(180) 


H = 


n x 0 Tly 

0 Tty n x 


(181) 


We let it be a finite element approximation to u. Furthermore, let G(u) denote a recovered 
approximation to be the stress obtained using any of the techniques we have just discussed. 

The so-called Cantin-Loubignac iteration is a process for further enhancing the accuracy 
of Cr(it) ~ a. The variational form of (179) is to find u such that 


where 


u E X : B(u,v) = (f,v) + g ■ v ds V v E X 

B(u,v) = J (SvYD(Su)dx 

if,v) = / f -vdx 

Jn 


and 


X = ju € [i/ 1 (fl)] 2 :u = 0 on r^j 
If our approximation to the stress is exact, then we would have 

J ( SvYG(u ) dx = J f ■ v dx + Jr g ■ v ds V v E X 
Thus, the lack of accuracy in G(u) is measured by the quantity 

(r,v)= f f-vdx+ f gv,ds— f (5v)‘G(w) dx 
Jti Jr w 

Since G(u) is continuous (by construction) we may integrate (187) by parts to obtain 


(182) 

(183) 

(184) 

(185) 

(186) 
(187) 
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(r , v ) = J {/ - S‘G(«)} -vdx + ^ {g~ HG(u)} ■ v ds V v G X (188) 

Notice that every term on the right hand side is computable. Furthermore, we have from 
(187) 

(r,v) = B(u,v) - [ (Sv)’G(u)dx = [ (Sv)‘{<r-G(u)}dx (189) 

v JQ JQ 

The method is therefore to first solve the approximate correction problem 

wex B(w,v) = (r,v) VveX (190) 

by means of a finite element discretization precisely the same as that used to obtain u. 
Having calculated w , we then set 


u* = u + w€ X (191) 

and compute <r* = G («*) in the same way as before. It is clear that <7* « cr. In fact, it is 
found that <r* is a much better approximation to the true stress than G(u). Of course, having 
obtained a*, the process can be repeated as often as we like. It is found computationally 
that the sequence of iterates tends toward <r, and that convergence of the iterates occurs in 
only a few iterations. 

The philosophy behind this approach appears similar to the Defect Correction Technique 
frequently used in numerical analysis. However, there are certain differences making the 
question of a mathematical analysis somewhat recalcitrant. In spite of this, the method 
works in practice. The major disadvantage associated with the method is the significant 
computational expense entailed while its major advantage is its generality. 

Extraction Methods 

A good overview of these techniques was presented in [39] which is summarized in this 
section. The concept of extraction formulas was first developed by Babuska and Miller 
[40,41,42]. They noticed that some properties of solutions u to a given boundary value prob- 
lem can be characterized by integral identities depending on u which they called extraction 
formulas, which are similar to Green’s formulas for elliptic problems. They proposed to cal- 
culate these properties with extraction formulas but with the exact solution u replaced by its 
FE approximation u. It turns out that with this method we can “extract” at a given point 
the value of u, derivative of u(x), its average value, etc., and the results are more accurate 
than those obtained by direct computation. 

Some spectacular examples showing improvement of accuracy due to applying extraction 
formulas were presented in [40] and [41]. Among them, the evaluation of stress intensity 
factors seems to be the most interesting (though, in this case, their approach is very closely 
related to other commonly known methods of linear fracture mechanics, for instance, [43]). 
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However, in none of the work in [41] was the problem of estimating gradients for two- 
dimensional problems considered at all, except for the calculation of normal derivatives on 
the boundary dfl of the domain of u. 

Here we follow the general idea of the method described in [40] and we develop appropriate 
formulas for extracting du/dx. We focus our attention on the model elliptic problem: 


d_ 

dxi 



u 


/(*) in ft 
u 0 on dfl 


(192) 


Integrating twice by parts the product of this equation and any function (tp + <p), we obtain 
the following identity: 


$ = J f • (y? + <p) dx -|- [- (ay?,,) , u ] + J Q [( a - a °) ^),«] u ^ x 


where 


$ 

$ 

$ 

a 0 

<P 


0 if (y? + (p) is smooth, at least a C 1 function, 
u(y) if <p{x, y ) = 7^ log \x = y[, 


2ira n 


du 

dx\ 


(y) if y ?(as,y) = 


1 


COS Q 


2na 0 I* - y|’ 


cos a = 


Xi - V\ 

I* - y I 


a(y), 

any smooth function (C 1 ) satisfying: 


(y? + <p) = 0 on dfl 

£-(<? + &) = 0 on dfl. 

on 

It is easily seen that in the second and third cases y?(*,y) was chosen as a solution of the 
given problem with the right hand side as Dirac delta at y and its generalized derivative, i.e., 
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dipole (but for a = a 0 = constant). That is why we recover u(y) and du/dx(y). The function 
<p was picked such that {ip + $) contains the appropriate singularity of but vanishes with 
its normal derivative on 5ft so that no boundary terms result from integrating by parts. The 
last term in our expression takes into account possible dependence of the coefficient a(x) 
on x and because of its strong singularity ( r ~ 3 ) it should be understood in the sense of the 
Cauchy Principal Value. 

The predicted rate of convergence of extracted quantities is equal to the double rate of 
convergence of u in the energy norm. This means that, for instance, for bilinear elements we 
should get ( h 2 ) error which is superior to that of directly computed derivatives. 

The method gives very good accuracy and has several attractive features. However, it 
has some important drawbacks that essentially limit its applications. 






it requires calculation of strongly singular integrals over the entire domain, 


there are difficulties with constructing a function (p for complicated domains; moreover, 
the more irregular the (p, the worse the quality of results given by the extractions 
procedure, 

there is a possibility to avoid the above two problems by restricting calculation to 
subdomains of a chosen regular shape (for instance squares); however, it turns out 
that such a procedure deteriorates the results very significantly; 


accuracy deteriorates also very rapidly as one approaches the boundary of ft, and 


a /'U in a/4 

wuvn 10 liiiuuou 


f ^ alliniip onor n f Avc fnr nirViir’ln 
uo viu puiv wpC'iMivv/iy '* mv*i 


know fund^mcntct! solutions 


(at least after averaging coefficients). 


Recovery of Boundary Fluxes 

As we have already mentioned, it can be very awkward to obtain good recovered approx- 
imations to fluxes near the boundary. The following method [44] is tailored toward deriving 
approximations to fluxes at the boundary. 

Consider the problem 


Suppose u € X C X = 
Thus 


-A u + u = f in 
u = 0 on 

{u € H'itl): v = 0 on 5ft}, is 



(194) 


the finite element approximation to u. 


u 


£X: J ( Vii • Vr + uv) dx = J fvdx 


VveX 


(195) 
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Let X° denote the set of shape functions which 
Then for any <pj £ A -0 we obtain from (194): 

/ ( Vu • Vtpj + wpj) dx = 

J n 

or rearranging we obtain 

/ = / (Vu • + uv?j) dx - [ fvjdx (197) 

JdQ an Jn Jn 

Replacing u by u on the right hand side gives the following statement for our approximation 
g ss on dCl: 


vanish at all the interior nodes on the mesh. 

/ f<pjdx+<£ Vj-z-ds (196) 

Jq Jd n on 


l <p jg ds= I (Vu • Vtpj + dx - I fa dx V^€X° (198) 

Jan J n Jn 

We now expand ^ in terms of the elements of A r °: g(s) = 2j = x ctk<pk{s) and substitute into 
(198) to obtain the system: 


Mot = b (199) 

where M is the n x n matrix with entries 

Mij = <f <pi{s)<pj(s)ds ( 200 ) 

JdQ 

and 6 is the n vector with entries 

6, = f (Vu- V<pi + iupi) dx - [ ftpi dx (201) 

Jn Jn 

In particular, 6, are rather easy to compute using the equations for the finite element ap- 
proximation u before the boundary conditions are imposed. It can be shown [45] that the 
resulting scheme is superconvergent. The scheme is very inexpensive owing to the presence 
of the tridiagonal matrix M. 

Summary and Conclusions 

We have described several methods for postprocessing. The remarks made in the text 
show that there is a great danger that several of these methods will not work for general 
approximations. The most promising methods seem to be those based on superconvergence 
and projection methods in conjunction with the procedure described Cantin-Loubignac it- 
eration. The use of the Cantin-Loubignac procedure is critical in the case of non-uniform 
meshes since otherwise the superconvergence techniques will fail. Meanwhile, using only 
projection techniques would result in very poor gradients near to the boundary. 

The main drawback of these is the cost, but it is conceivable that by working locally one 
could significantly reduce the costs without sacrificing accuracy. 
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x — 0 £2 j ^2 £2 ^ ♦ . . . x = 1 
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► 


h 


Nh =1 


Figure 3.44: Uniform mesh in R 1 . 

Extraction methods are rather unsuited for use with a general purpose package owing 
to their high dependence on the features of each particular problem. Moreover, they are 
expensive to implement and to use. 

Boundary flux recovery methods represent a very attractive method for dealing with 
recovery of fluxes on the boundary and again is valid with general types of meshes. Thus the 
most favorable techniques can be classified according to the portion of the domain in which 
they are to be used: 

1. For postprocessing on the interior of the mesh, use superconvergence and projection 
methods in conjunction with Cantin-Loubignac iteration and experiment to see if they 
can be applied locally. 

2. For postprocessing on the boundary, use the boundary flux recovery technique. 
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Figure 3.45: Postprocessing when p = 1. 
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Pass 


Fail 


Figure 3.47: Meshes satisfying and failing the “six triangle” property. 
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Figure 3.49: Regular mesh of quadrilaterals. 
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(a) 


PROJECT: Pure Bending DEFORMATION - 5x Magnification 



Figure 3.50: Pure bending of a prismatic beam, (a) Beam with applied moment, (b) deformed 
shape of beam (magnified five times). 
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PROJECT: Pure Bending Zienkiewia-Zhu Error Eann«tt 



ERRORS 
MIN=35.003712 
MAX=65 .48594* 
GLOBAL=82.09098 1 


(a) 



PROJECT: Pure Bending Amiytic Eiror 


ERRORS 
MIN-66.433897 
MAX=78.00471 
GLOBAL* 128.77473 


(b) 


Figure 3.51: Pure bending of a prismatic beam analyzed using a mesh of three linear ele- 
ments. (a) Error estimates calculated using the Zienkiewicz-Zhu method, (b) analytic errors 

(measured in the energy norm). 


PROJECT: Pure Bending Zienlriewici-Zhn Error Efiimaie 



ERRORS 
MIN* 11.925241 
MAX*14.813858 
GLOBAL^ .01 8822 


(a) 



PROJECT: Pure Bendmj Amiytic Error 


ERRORS 
MIN=12.1*4312 
MAX« 14.237051 
GLOBAL=66-55Q5 16 


(b) 


Figure 3.52: Pure bending of a prismatic beam analyzed using a mesh of 24 linear ele- 
ments. (a) Error estimates calculated using the Zienkiewicz-Zhu method, (b) analytic errors 
(measured in the energy norm). 
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PROJECT: Pure Beading 


Zienkiewicz-Zbu Error Fgrm ai r 



ERRORS 
MBW.8929546 
MAX«2.6963 307 
GLOB AL=3 1 .37 1 13 8 


(a) 


PROJECT: Pure Bending 


Analytic Error 


ERRORS 
MIN=2.0807615 
MAX=2J5 067306 
GLOB AL*3 2.47 1543 



(b) 


Figure 3.53: Pure bending of a prismatic beam analyzed using a mesh of 192 linear ele- 
ments. (a) Error estimates calculated using the Zienkiewicz-Zhu method, (b) analytic errors 
(measured in the energy norm). 
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PROJECT: Parc Bending 


Zeakiewicz -Zha Error E -s nm aie 



ERRORS 
MIN *0323 7849 
MAX-0.482018 
GLOBAL- 15.785832 


(a) 


PROJECT: Pure Bending 


Analytic Error 


ERRORS 
MIN-03668845 
MAX-0.4434079 
GLOBAL- 16. 10589 



(b) 


Figure 3.54: Pure bending of a prismatic beam analyzed using a mesh of 1536 linear ele- 
ments. (a) Error estimates calculated using the Zienkiewicz-Zhu method, (b) analytic errors 
(measured in the energy norm). 
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3.3.2 Error Estimate Application: The Zienkiewicz-Zhu Technique 

While the methods in the previous section were examined in the context of solution enhance- 
ment, they can also be applied to error estimation. Namely, an enhanced solution can be 
compared to the original solution to provide an estimate of the error in the original solution. 
For example, the Zienkiewicz-Zhu family of error estimates is based on superconvergence 
and projection methods. In this section, we present the results of a study of an enhanced 
version of the simple error estimator proposed by Zienkiewicz and Zhu [16] and analyzed by 
Ainsworth, et al. [38]. This technique was incorporated into an existing /ip- adaptive finite 
element analysis code for testing. Following an outline of the technique, a few examples are 
shown for several types smooth problems. 

For a governing differential equation and boundary conditions of the form 


—V • (AVu) + B ■ Vtt + Cu = / in ff 

u = /o on da 

an estimate of the element-wise error is sought in the following semi-norm: 


E (f Ik) = ^ (Vu — Vit/,) • [A (Vu - Vu*)] 


d ft 




The error estimates are calculated by first calculating the finite element solution, u* G V h , 
and then projecting the derivatives of u* onto V* to obtain G*. Namely, Gh G Vh x Vh satisfies 

irun f [A(F h - Vu k )] • (F h - Vu h )dQ 

The Zienkiewicz-Zhu error estimate is then calculated according to: 

l 

E h (Sl) = fc / - Vu a) • [A (G h - Vti*)] dQ 2 

L *: J{1 * J 

Sample Results 

In order to compare the error estimates to the analytic errors, the effectivity indices are 
defined globally and elementwise as 

_E h (a) _E h (a k ) 

7n E(a) E(a k ) 


where E(Q) = 




1 

2 
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Example: Problem With Polynomial Solution 

For this example, the loading function / and the boundary condition f 0 were chosen 
to correspond to the simple polynomial: 

+(»-!) “I 

Figure 3.55 shows a contour and a carpet plot of this function along with the finite 
element mesh used in this example. Figures 3.56, 3.57, and 3.58 show plots of the 
elementwise effectivity indices and the estimated error for the following forms of the 
differential equation: 


Figure 3.56: 


_V . Vu = / 


Figure 3.57: —V • Vu + u = / 

_ du 

Figure 3.58: -V ' Vu + ^ = / 

As shown in these three figures, the estimated errors are quite accurate (the effectivity 
indices are close to one) for all three forms of the governing differential equation. 

Example: Problem With Sinusoidal Solution 

ih ^ | V, n fun/'t irm f onrl f Vio Kaii n rjan? COndlvlOU- fcii W6F6 chosen 

lGr l ni s example, uut lvyu-um^ iuuvwvu, j , ^ *****■' j JU , • 

to correspond to the analytic solution: 


u = sin(x + 3y) + e x + y 2 

Figure 3.59 shows a contour and a carpet plot of this function along with the finite ele- 
ment mesh used in this example. Figure 3.60 shows plots of the elementwise effectivity 
indices and the estimated error for the differential equation: 


—V • Vtt = / 

Except for a few elements, the estimated errors in this example are also quite good. 

Example: Problem With Arctangent Solution 

For this example, the loading function, /, and the boundary condition, / 0 , were chosen 
to correspond to the analytic solution: 

u = tan -1 (15x) 
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Figure 3.61 shows a contour and a carpet plot of this function along with the finite 
element mesh used in this example. In Fig. 3.62, plots of the elementwise effectivity 
indices and the estimated error are shown for the differential equation: 

-V • Vu = / 

For this example, the estimated errors were not as accurate as the previous two exam- 
ples, possibly due to the limitations of the technique. 

Example: Elasticity problem 

For this example, a short cantilever beam was analyzed using the sequence of meshes 
provided in Ainsworth, et al. [38]. Figure 3.63 shows the geometry and loading of 
the cantilever beam while Figs. 3.64 and 3.65 show the sequence of six /i-adapted 
meshes. These meshes were used to illustrate another (apparent) feature of this type 
of error estimate. Namely, it is possible to estimate the effectivity index of a given 
error estimate. This is done by first computing the finite element solution and error 
estimate on the original mesh and then globally h-refining the mesh and recalculating 
only the error estimate on the newly-refined mesh. The effectivity index appears to 
obey the relationship: 

7 ( 

where r « 3. 

Furthermore, this relationship can be used to produce an improved error estimate 
according to: 


Eh (initial mesh)' 
E 


Eh (refined mesh) 


[ Eh (initial mesh) J 


r, ,. 1N [^(original mesh)] r+1 

£„(,mproved) = | £)t(refined mesh)r 

Figure 3.66 shows the effects of these relationships for the six meshes of the cantilevered 
beam problem. In this figure, the estimated errors and effectivity indices are shown 
for various values for the exponent r. Namely, curves “estim” and “effec0” correspond 
to the unmodified error estimates while “estnewl” and “effecl” correspond to a value 
r = 1, “estnew2” and “effec2” to a value r = 2 and “estnew3” and “effec3” to a value 
r = 3. As shown in this figure, this technique can greatly enhance the quality of the 
error estimate although the added expense of obtaining solution and error estimates 
on additional meshes makes the technique impractical. 
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DJECT: new 


SOLUTION 


PROJECT: new 


SOLUTION 



| i i ] ^ ;n 

^1 4S AS *2 - 0.3 


MIN=-7.749V93 
MAX=-0. 750001 



MTN^774V993 
MAX =-0.750001 


PROJECT: new MESII 



Figure 3.55: Plots of analytic solution, u = (x - *) 3 + (y - 5) - 1 and the finite element 
mesh used in the first example. 
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PROJECT: new 


LOCAL EFFECTTVITY INDEX 







, r ■ 

t 

•r " - 

1 1-5 

2 

PROJECT: new 


ERROR ESTIMATION 



0.2 °£ 



effectiviry 
MIN=0.6716723 
MAX= 1.23 1825 
GLOBAL=0.9616329 
operator 
-V(V u)=f 
D.O.F= 289 


MIN =0.0285 68 
MAX=1. 0245753 


Figure 3.56: Local effectivity indices and estimated error for example with polynomial solu- 
tion with the governing equation: 


-V • Vu = / 
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effectivity 
MIN=0.6722238 
MAX=1.2293515 
GLOBAL=0.9625260 
operator 
-V(Vu)+u=f 
D.O 289 


PROJECT: new 


ERROR ESTIMATION 



MIN=0.0285986 
MAX=1. 0254241 


Figure 3.57: Local effectivity indices and estimated error for example with polynomial solu- 
tion using the governing equation: 


— V • Vu -f u = / 
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PROJECT: new 


LOCAL EFFECUVTTY INDEX 



effectivicy 
MIN=0.6724334 
MAX-1.2332381 
GLOBAL=0.96 11406 
operator 
-V(Vu}+du/dx=r 
D.O.F= 289 


PROJECT: new ERROR ESTIMATION 




MIN=0.0285454 
MAX=1. 0259369 


Figure 3.58: Local effectivity indices and estimated error for example with polynomial solu- 
tion using the governing equation: 


_ _ du 
-V-Vu + ^-Z 
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PROJECT: new 
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Figure 3.59: Plots of the analytic solution, u = sin(x + 3 y) + e x + y 2 and the finite element 
mesh used in the second example. 
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PROJECT: new 


LOCAL EFFECTIV1TY INDEX 



effectiviry 
MIN=0.2792334 
MAX-3.2S51514 
GLOBAL=0.9893719 
operator 
-V(Vu)=f 
D.O J= 289 


PROJECT: new ERROR ESTIMATION 



MIN- -0.453e-03 
MAX-0.6711019 


Figure 3.60: Local effectivity indices and estimated error for the second example with the 
governing equation: 


-V • Vu = / 
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PROJECT: new MESII 
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Figure 3.61: Plots of the analytic solution, u = tan ^lSx) and the finite element mesh used 
in the third example. 
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PROJECT: new 


LOCAL EFFECTTVITY INDEX 



effecuvicy 
MIN=0. 12026 13 
MAX=2.6979657 
GLOBAL=J44882 
operator 
-V(VuW 
D.O Ja 289 


PROJECT: new 


FIRST COMPONENT 


* - 0.1 




0.03 0.25 03 0.63 


MIN=-0.053911 

MAX=0.6154387 



Figure 3.62: Local effectivity indices and estimated error for the third example with the 
governing equation: 


-V • Vu = / 
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PROJECT; c«nl MESH 6 



Figure 3.65: Meshes 4, 5, and 6 for the cantilever beam problem. 
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Error estimation for four estimates 



Effectivities of four estimates 


Figure 3.66: Error estimates and effectivity indices for the 6 meshes for the cantilevered 
beam example. Shown are results using different exponents for r for improving the quality 
of the error estimates. 
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3.3.3 Error Estimation and Solution Enhancement: The Projection Method 
with Cantin-Loubignac Iteration 

In the previous section, error estimation results were presented using the technique proposed 
by Zienkiewicz and Zhu [16] and analyzed by Ainsworth, et al. [38]. This technique was 
used to apply projection methods to Poisson problems and problems in elasticity in order 
to obtain improved stresses as well as estimates of the error in the approximate solution. 
In addition, Section 3.3.1 introduced a technique called Cantin-Loubignac iteration which 
can be used to improve the basic error estimate. This technique is discussed further in this 
section. 

First we recall that we seek a solution to the differential equation 

- V • (AVu) + B ■ Vu + Cu = f (202) 

in the domain ft, subject to the boundary conditions 


= ii 0 on T d 
= t on Tt 

such that T/j U F( = 5ft and n is the outward-directed boundary normal. 

The variational form of this problem may be stated: Find u E Vo(ft) + {«o} such that 

B{u,v) = F(v) VueFo(ft) (203) 

where Vo(ft) = {v: v € if 1 (ft), *yv = 0 on T^} and 


u 



B(u,v) 

F(v) 


f [AVu • Vu + B ■ Vuu + Cuv]dCl 
Jn 

f /udft-f / tvdS 
J n Jr, 


(204) 


The basic Zienkiewicz-Zhu error estimate is obtained by first solving the approximate prob- 
lem for the finite element solution, Uh 


The projection problem 


B(u h ,v h ) = F(v h ) 
is then solved to obtain 


Vi>* 6 V k C V 0 

the enhanced solution derivatives, 


(205) 


G h 




Vv h e v h x v. 


(206) 
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and the Zienkiewicz-Zhu estimate is obtained 


e „{ n) = 



(G k -Vu h )-A(G h -Vu h )dQ 


1 

2 


(207) 


This estimate can be improved by using the enhanced derivatives to construct a correction 
problem 


B (6, v h ) = F K) - / (AG h • Vv fc ) (208) 

J Q 

for 8 € Vh and V Vh € V&. 

The approximate solution is then enhanced by adding the correction 

ut =u h + 6 (209) 

If the size of the correction, ||6||, is not small, the enhanced solution derivatives, G&, are 
recalculated using equation (206) with the updated solution, uj, and improved estimate 
using equation (207), and then a new correction using equation (208). This cycle can be 
repeated until ||6|| is sufficiently small. 

As an example of this technique, we analyzed the plane elasticity problem of a short can- 
tilevered beam from the previous section using the sequence of meshes provided in Ainsworth, 
et al. [38]. Figure 3.63 shows the geometry and loading of the beam while Figs. 3.64 and 3.65 
show the sequence of six /^-adapted meshes. The above iterative procedure was performed 
until ||6|| was six orders of magnitude smaller than ||tt + ||. Typically, only six or seven iter- 
ations were necessary. Figure 3.67a compares the modified error estimates to the standard 
Zienkiewicz estimates and actual errors (taken from [38]) for the six meshes. Similarly, Fig. 
3.67b compares the effectivity indices (error estimate/actual error) for the two estimates. 

As these figures show, the modified technique provides an improved estimate in all of the 
cases. Unfortunately, each iteration of this technique requires the solution of a correction 
problem that is just as expensive as the original solution. Thus, with the added cost of the 
projection problems, this technique is impractical on a global (meshwise) scale. In the next 
section, we examine the consequences of using this technique on a patchwise basis. 
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3.3.4 Localization Study of Solution Enhancement Techniques 
Introduction 

The purpose of this work was to examine the various types of postprocessing schemes 
suggested by a survey described in Section 3.3.1. In particular, the suggested techniques 
were studied in the context of finite element approximations and applied to one element at 
a time instead of to the entire mesh. 

In order to simplify our explanation, we employ a very simple problem to illustrate the 
techniques. 

Let 0 C E 2 be an open bound domain with Lipschitzian boundary 5fL For given data 
/ € L 2 ( 0), and g € L 2 (r^) we seek the solution it of 


— V 2 u = / in ft 

u = 0 on Tp 
du 

^ = sonr " 

where T& fl I"V is empty, Tp is non-empty and dfi = Td U Tw- 
The variational form of (210) is to find 


( 210 ) 


u€ X : B(u,v) = L(v) VveX (211) 

where 


X = {u € if^O) : 7U = 0 on T/)} 
B(u,v) = fn'Vu-'Vvdx > 

L(v) = fa fvdx + f Vs gvds 


( 212 ) 


Let X C X denote a space of finite elements defined on a partitioning P of 0 into the union 
of subdomains {0 a }^ =1 satisfying the usual assumptions. The finite element approximation 
u « u is defined by 


u € X : B(u , v) = L(v) Wv e X (213) 

Having computed u, it is usual to postprocess u to obtain some approximation to the gradient 
which we shall denote by p « Vu. 

In the following sections, we consider various schemes by which p can be obtained. For 
example, one may perform some projection of Vu onto smooth functions (see Hinton and 
Campbell [35], Rachowicz and Oden [36]), or, at the opposite extreme, simply take p V«. 

The study described in these sections attempts to answer the following questions. 
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Short Cantilever Beam : Error Estimates 



Short Cantilever Beam : Effectivity Indices 



(b) 

Figure 3.67: (a) Exact error and two error estimates for the six meshes of the cantilever 
beam problem, (b) Effectivity indices of the two estimates. 
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• how accurate is the approximation p « Vu and can we estimate the accuracy numer- 
ically? 

• how accurate is the approximation p ss Vu locally and can we perform some local 
procedure by which to enhance the local approximation properties? 

The first of these questions really addresses a much wider issue than is immediately 
apparent. Namely, the problem of designing a general purpose error estimation package and 
solution enhancement procedure. 

The basic purpose of such a package is to accept some approximation p to the true gradi- 
ent of the solution of a given boundary value problem and to produce a numerical estimate 
of the accuracy p « Vu and, possibly, an improved approximation over p. Therefore, the 
first question is essentially how to design an error estimation and postprocessing package. 

The second question recognizes that while many postprocessing schemes appear satis- 
factory when viewed at a global level, they may, nevertheless, perform poorly at a local 
level especially near boundaries of the domain or near singularities in the true solution. The 
question then arises of whether and how the approximation p « Vu may be enhanced in 
these areas. 

A General Framework for Solution Enhancement and Accuracy Assessment 

In this section we outline how one can mathematically formulate the problem discussed 
in the introduction. Suppose, once again, that p « Vu is some postprocessed gradient. Let 
the error in this approximation be denoted by E. That is, 

E = Vu — p (214) 

Obviously, we cannot compute E exactly. However, it is worthwhile to note that if we could 
compute E (or at least a good approximation to E) then we not only have a means of 
computing the accuracy of p « Vu, but also a means of improving p, thereby obtaining an 
enhanced approximation. That is, if E* ~ E then 

Vu =£ + £«£ + £* (215) 

and, consequently, p + E* is an enhanced approximation. 

The question now, is how can E be approximated since Vu is unknown? 

Let v 6 X (where X is given by (212)). Then, multiplying (214) by Vu and integrating 
over fi gives 


X 


E • Vu dx 


/ Vu • Vu dx 

Jn 


[ p • Vu dx 

Jn 


Now, recalling (211) and (212), we obtain 


( 216 ) 
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/ E ■ Vu dx 
J n 


B(u,v) — p • Vu dx 

J o 

L(v) -Ip- V« dx 

Jn 

J fv dx + J gv ds — 



■ ' Vv dx 


That is, for any v € X 

f E • Vt> dx = f fv dx + [ gv d s — f p ■ Vw dx (217) 

J fi Jq Jr n Jri 

Examining (217), it is clear that every term on the right hand side can be computed. There- 
fore, it is possible to calculate E using (217) and hence obtain accuracy estimates and 
enhanced approximation as described above. While this approach is clear and straightfor- 
ward, it does suffer an important drawback: it requires the approximation and the solution 
of a global problem and is therefore expensive. 

One issue which we shall investigate is that of how to decompose (217) into a sequence 
of local problems on each element in the partition P. Solving a sequence of local problems 
is significantly less expensive than solving a single global problem. 

Analysis of Nodal Averaging Postprocessing Schemes 

We now take up the problem described in the previous section for the specific case of 
p being a simple nodal average of gradients. For simplicity, we assume Q (the domain) is 
rectangular and has been partitioned into the union of squares of side h. The space X is 
taken to consist of piecewise bilinear finite elements. Figure 3.68 shows how p is defined at 
the nodes in the partition: for an interior node x a, we define 


p{x A ) = l - { Vu|cj + Vu|c 2 + Vu|c 3 + Vu| c< } (218) 

where Vu|ci denotes the direct approximation to the gradient evaluated at the centroid Cj 
of element Qj. For a boundary node *jwe define 


p(*fl)= 5 (Vfi|c s +Vi|c.} (219) 

where and fl 6 are shown in Fig. 3.68. 

Proceeding in this manner allows us to obtain postprocessed values at each node in the 
mesh. Having done this, the value of p at other points is obtained by means of bilinear 
interpolation of the nodal values. 

This method is not being proposed as a postprocessing technique although it does coincide 
with a recognized superconvergent scheme on the interior of Q, (see Ainsworth and Craig 
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[15], Zlamal [46]). Since it represents a sensible choice for p, however, it has been chosen to 
illustrate our ideas. Furthermore, since this scheme is inferior (non-superconvergent) near 
the boundary, it represents a good choice for testing whether our methods for postprocessing 
can correctly identify the poor approximate gradients near the boundary and correct them 
accordingly. 

Example: 1 

As a first test problem we consider 


— Vu 

u 

du 

dn 


0 in fl 


x 3 — Zxy 2 on 

/ 3x 2 - 3 y 2 \ 
H ' [ - 6 xy ) 


on Tw 


where fl,rAr and T/j are shown in Figure 3.69. The true solution to the problem is u = 
x 3 — 3 xy 2 , meaning that we can compare our numerical approximations with the true solution. 

In order to compare our subsequent local schemes with the accuracy attainable without 
localization, we first solve (217) globally, using the so-called Cantin-Loubignac Iteration, 
which we now briefly describe. First, we replace E by V<^» in (217) (where (j) G A), giving 


(j> G X : / V<f>- Vt> dx = / fv dx + / gvds — / p • Vu dx Vu G X 
Jq Jq Jr N Jq 


( 220 ) 


where p is the gradient obtained as described above. Integrating by parts allows us to write 
(220) in the form presented by Cantin and Loubignac originally: 


(f> G X : / V<f> ■ Vv dx = / (/ + V • p)v dx -f / ( g — n ■ p)vds Vu G X 
Jq Jq Jr ^ 


(221) 


We now discretize (221) in the same way as the original problem was discretized: 
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4> € X : / V<£ • Vu </x = / (/ + V • p)v dx + / (y - n • p)u ds Vu € X (222) 
Jn do -/rA, 


This problem is then solved in the same way as was u. Having obtained we set E = V<^ « 
£ and perform the correction 


u new = u + <£ 
p new = G(u new ) 


(223) 


where G(u new ) denotes the gradient approximation obtained by performing nodal averages 
of u new . Naturally, we can now use p new in place of p in (222) and repeat the procedure. 
Table 4 shows the accuracy of these iterates along with the accuracy of Vu ~ Vu and 
p fa V« (i.e. when no correction has been performed). The following features can be seen: 


• the uncorrected approximation p « Vu is inferior to not applying any smoothing. 
That is 

||Vu - Vu|| L 2 (n) < ||Vu - p|| L 2 ( fi ) 

The reason for this can be seen by examining Table 5 which shows the contributions 
to the global error from each element. The effect of the nodal averaging is to degrade 
the approximation in elements adjacent to the boundary. 

• the sequence of Cantin-Loubignac iterates provided a sequence of improved approxima- 
tions to Vu. These iterates converge in only a few steps and the greatest improvement 
is seen in the first application of the procedure. 

• examining the local behavior of the iterates (Table 5), we see that the method correctly 
identifies the poor approximation of p along the boundaries and improves the accuracy 
there. 


In summary, the Cantin-Loubignac scheme works successfully, although the actual results 
are degraded owing to the deficiencies in the nodal recovery scheme being employed. Finally, 
we might ask what is the rate of convergence of the Cantin-Loubignac iterates as the partition 
is refined? Table 6 shows the results obtained on progressively finer meshes, from which it 
appears that the rate of convergence is slightly better than <r(h). 

Localization Strategies 

Let us now return to our earlier question, namely that of how to localize the problem 
(217). Proceeding from (220) we have 
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(224) 


4>€X : J ^7(j) ■ Vu dx = J fv dx + J gvds — J p • Vu dx Vu € X 


To directly localize (224) will result in a very poor scheme. Therefore, we first rewrite (224) 
by integrating by parts and then decomposing into a sequence of problems over each element 
in the partition: for any u 6 X, 


L fvdx+ L gvds — J p ■ Vu dx 


f (/ + V -p)v dx + f (g — n ■ p)vds 
J q J r N 

(g — n ■ p)vds 


/an^-nr* 


E {/„(/ + v-p)w* + / a 

V / fvdx+ gvds — p ■ Vu dx 

^{Jqk Jda K n r * Ja K 


) 


+ / (ni< ■ p)vds \ 

Jdn K \v N J 

In view of (225), we decompose (224) into the sequence of local problems: 


fc\c\r \ 


4>k € X K : / V<t> K -Vvdx = f n fv dx + gvds 

J Q h - JdQ K C\ r N 

- fn K P • Vvdx + fan K \ r N n K • pvds (226) 

Vu e x K 

where Xk = {u € H l {Slx) : ■yv = 0 on dVt K DTi)} and n h denotes the unit outward normal 
vector on dflh - 

One advantage of the localization process is that we can calculate a high order approxima- 
tion (j>K to the solution of this problem. That is, we need not restrict ourselves to calculating 
a bilinear approximation. However, there are drawbacks to the localization. In particular, 
the boundary conditions for the local problem are imprecise, whereas they are known ex- 
actly for the global problem. Therefore, the major difficulty in localizing the problem is the 
determination of the boundary conditions on edges which lie on the interior of the mesh. 

Example: 2 
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The purpose of this example is to investigate the influence of various choices of boundary 
conditions on the interior edges. Once again we use the simple test problem 

— Vu = 0 in 
u = u on 

du 

d ^ =9 ° n 

where I'd, Tn are shown in Figure 3.69. As before the true solution is u(x , y) = x 3 — 3 xy 2 in 
n. In Table 7 we compare the accuracy obtained by solving the local problems with three 
choices of boundary conditions. A typical local problem stated in strong form is 

Find <f>K- 



(227) 


— V (f>K 
d<t>K 

dn 

<f>K 

d(f>K 

dn K 

where g * should be chosen as an 
for the following choices: 


= / + V-p 

= 9 -n p 
= u — u 

= 9’ ~n K p 


in Qk 

on 

on 

on dClh \T 


(228) 


approximation to the true flux fix • Vu. We show results 


(a) 

* du 

9 — dn K 

(i.e. using true flux) 

(b) 

9* = n K 'P 

(i.e. using approximation to true flux based 
on nodal averaging) 

(c) 

9’ = \n K ■ (Vu|i eft + Vfi| right) = 

(t^) (i.e. the flux approximation obtained by av- 


eraging fluxes between neighboring elements) 


The first case is not a practical method, but is included to show the best results which 
can be attained. Examination of Table 7 shows that choice (b) performs rather poorly due 
to the behavior in elements along the Dirichlet boundary. Better results are obtained using 
choice (c) which is almost as good as using the true fluxes 

A more significant conclusion is reached by comparing the accuracy of the local scheme 
with the full Cantin-Loubignac iteration. Referring to Table 4, one sees that at convergence 
the accuracy of the full Cantin-Loubignac scheme is .298 compared with .201 obtained using 
the local scheme. The cost of the local scheme is considerably less than the full scheme and 
is more accurate. 

The conclusions which we may draw from this simple investigation are: 
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1. not only is a localized version of a Cantin-Loubignac scheme viable, but in fact, per- 
forms better than the global scheme (owing to the possibility of solving local problems 
using high order methods). 

2. the choice of boundary condition for the local problems is critical. 

To enlarge upon the second point, we have seen that neither the Cantin-Loubignac scheme 
or the localized version can fully compensate for the deficiencies in a poor choice of post- 
processing scheme to obtain p. However, they can locate the areas of the mesh where the 
scheme is poor. 

Further Numerical Studies 

In this section we continue the study outlined previously. However, we now consider an 
alternative choice of postprocessing by which p can be obtained. Previously, we employed a 
deliberately poor scheme. 

The choice of p is similar to before, except that we now do not choose centroid values 
(c.f 218-219) but choose to average values at the node. That is (c.f. Figure 3.68) 

p(*a) = j {Villain, -i- Vu| Xy4 ,n 2 + Vu| Xy4 ,n 3 + Vu| Xyt ,n 4 } (229) 

where ^u\x A ,Ui denotes the value of Vu obtained by evaluating Vu on element at 
global coordinate x^. The values at other nodes are determined in a similar way. Values at 
other points are obtained by bilinear interpolation of the nodal values. This scheme is not 
superconvergent, but once again represents a sensible choice of postprocessing. 

Example: 3: Full Cantin-Loubignac Scheme 

Once again we solve the simple problem of the previous examples, but this time we use 
the full Cantin-Loubignac scheme in conjunction with the recovery scheme described above. 
The accuracy of the Cantin-Loubignac iterates are shown in Tables 8 and 9. 

Comparing these results with Tables 4-5, it is seen that the new choice of p leads to a 
significant improvement in the accuracy of the postprocessed gradient and the iterates. 
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Figure 3.69: Domain, mesh and boundary conditions for the first example 


Table 4: 
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.303001 
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.300554 
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.298661 


||e||£; = .498827 


Accuracy of Cantin-Loubignac iterates for the first example. 
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Iterate 0 





rror in 
Row: 4 

1 : o.: 

0.425E-02 

=2406025288 

0.392E-02 

887 

0 . S31E-02 

0 . 696E-02 

0 . 727E-02 

0.596E-02 

0.100E-01 

0.231E-01 

Row : 3 

0 . 335E-02 

0.298E-02 

0 . 158E-02 

0 . 55 6Z-03 

0 . 459E-03 

0 . 994E-03 

0.737E-03 

0.531E-02 

Row : 2 

0.179E-02 

0 . 832E-03 

0 . 632E-03 

0.343E-03 

0 . 383E-03 

0 . 803E-03 

0 . 484E-Q3 

0.519E-02 

Row: 1 

0 . 167E-02 

0 . 999E-03 

0 . 671E-03 

0 . 437E-03 

0.532E-03 

0 . 113E-02 

0 . 929E-03 

0.345E-02 


Table 5: Element by element errors for iterates in Cantin Loubignac Scheme for the first 
example. 
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Iterate 

1 Vu -p|U 2 (n) 

0 

.186232 

1 

.132681 

2 

.124869 

3 

.121490 

4 

.119525 

5 

.118210 

6 

.117248 


.116498 


.115886 


Iterate 

WIWJIWI 

0 

.066187 

1 

.048201 

2 

.045755 

3 

.044704 

4 

.044098 

5 

.043694 

6 

.043401 


.043176 


.042997 


v (i — v 


u llL 2 (n) 

h=i 


a An w r% 

= .^yooo 


ii vu- vuji L2(n) = 



.124982 


Table 6: Accuracy of Cantin-Loubignac iterates for the first example on finer meshes. 
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t»l t’l 


ERRORS IN DIRECT GRADIENT APPROXIMATION 
Error in H~l: 0.49882675632465 . 


Row: 

4 

0 . 180E-01 

0 . 924E-02 

0.338E-02 

0 . 452E-03 

0 . 452E-03 

0.338E-02 

0; 924E-02 

0.180E-0l| 

m 



0 . 924E-02 

0.338E-02 

0 . 452E-03 

0 . 452E-03 

0.338E-02 

0.924E-02 

QBS5S9 

Row: 


0 . 180E-01 

0 . 924E-02 

0.338E-02 

0 . 452E-03 

0.452E-03 

0.338E-02 

0 . 924E-02 

0.18OE-01 

m 



0.924E-02 

0.338E-02 

0.452E-03 1 



[o. 924E-02 

0.180E-01 


ERRORS IN CORRECTED GRADIENT: 
Error in H~l: 0.17647157085140 


Row: 

4 

0 . 540E-05 

0.949E-03 

0 . 949E-03 

0 . 949E-03 

0 . 949E-03 

0.949E-03 

0 . 949E-03 

0.174E-01 . 


B 



0 . 911E-05 

0 . 911E-05 

0 . 911E-05 

0 . 911E-05 

0 . 911E-05 

0.366E-02 


B 



0 . 911E-05 

0.911E-05 

0 . 911E-05 

0 . 91 IE-0 5 

0 . 911E-05 

0.208E-02 

Row: 

B 


0.216E-04 

0.216E-04 

0 .216E-04 

0.216E-04 

0.216E-04 




(a) Results obtained with exact boundary conditions. 


ERRORS IN CORRECTED GRADIENT: 


rror in 


0 . 100E-01 

0.100E-01 

0 .100E-01 



0 . 100E-01 

0.314E-01 

— 

Hi 

0 . 917E-02 

0 . 732E-04 

0.732E-04 

0 .732E-04 

0 .732E-04 

0 . 732E-04 

0 . 732E-04 


Row: 2 

0 .514E-02 




0 . 732E-04 

0.732E-04 

0 . 732E-04 

0 . 317E-02 

m 


0 . 857E-04 

0 . 857E-04 

0 .857E-04 

0.957E-04 

0 . 857E-04 1 

0 . 857E-04 j 

0 . 177E-02 


(b) Results obtained using n ■ p as a boundary condition. 


RRORS IN CORRECTED GRADIENT: 

rror in H'l: 0.20186638019028 A. 


Row : 

4 

0. 678E-03 

0 . 120E-02 

0.120E-02 

0.120E-02 

0.120E-02 

0.120E-02*0.120E-02 

0 . 177E-01 

Row: 

3 



0.256E-03 

0.256E-03 

0.256E-03: 

0.256E-03 

0.256E-03 

0.390E-02 

Row: 

2 

0 . 849E-03 

0.256E-03 

0.256E-03 

0.256E-03 

0 .256E-03 j 

0.256E-03 

0.256E-03 

0 .232E-02 

n 

B 

0 . 678E-03 



0.269E-03 

0 .269E-03 j 

0.269E-03 

0.269E-03 

0 . 160E-02 


(c) Results obtained using {§£) as a boundary condition. 


Table 7: Impact of different choices of boundary condition for localized problems. 
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Iterate 

Vu-p z, 2 (n) 

0 

.373434 

1 

.220942 

2 

.190204 

3 

.174521 

A 

i n a p * p - r\ 

.lV^OOZ 

5 

.157789 

6 

.153111 

7 

.149872 

8 

.147641 

9 

.146119 


Table 8: Accuracy of Cantin-Loubignac iterates for the example problem. 
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fj 1*1 


BOORS IN DIRECT GRADIENT APPROXIMATION: 
RKUrw i* * aa o') cn e crn a cz 


-r?r is-. 

Row: 4 

H~.ll 2-J 

0.180E-01 

0 . 924E-02 

0.338E-02 

0 . 452E-03 

0 . 452E-03 

0.338E-02 

0 . 924E-02 

0.180E-01 

Row: 3 

0.180E-01 

0 . 924E-02 

0.338E-02 

0 . 452E-03 

0 . 452E-03 

0.338E-02 

0 . 924E-02 

0.180E-01 

Row: 2 

0.180E-01 

0 . 924E-02 

0.338E-02 

0 . 452E-03 

0.452E-03 

0.338E-02 

0 . 924E-02 

0.180E-01 

SOI 

0 . 180E-0I 

0 . 924E-02 

0 . 338E-02 

0 . 452E-03 

0 . 452E-03 

0 . 338E-02 

0.924E-02 

0.180E-01 


ERRORS IN NODALLY AVERAGED GRADIENT: Iterate 0 


:ror 

Row: 

in 

4 

H A 1: 0 .; 

0 . 169E-Q1 

3734344982c 

0.521E-02 

1838 

0 .228E-02 

0 . 818E-03 

0.818E-03 

0.228E-02 

0.521E-02 

0.169E-01 

Row: 

3 

0.790E-02 

0.574E-03 

0 . 574E-03 

0.S74E-03 

0 .574E-03 

0.S74E-03 

0.574E-03 

0.790E-02 

Row: 

2 

0.790E-02 

0.574E-03 

0 . 574E-03 

0.S74E-03 

0.574E-03 

0.574E-03 

0.574E-03 

0.790E-02 

Row: 

3 

0 . 169E-01 

0 . 521E-02 

0.228E-02 



0.228E-02 

0.521E-02: 

1 

0 . 169E-01 


"RRORS IN NODALLY AVERAGED GRADIENT: Iterate 1 

- - -*'"‘ 2713 

|0 . 424E-03 

0 .164E-03 


-sr m h~±: 

Row: 4J0.432S-02 

0 .147E-02 

Row: 3- 

0.256E-02 

0.129E-02 



0 . 129E-02 

.EBB 

0.432E-02 

0 .147E-02 


0 . 150E-03 jO . 148E-03 
0.737E-04|0.725E-04 


0 . 164E-03 


0.424E-03 0 


0 . 737E-04 lo.725E-04 |0.152E-03 
150E-03"|o .148E-03 


0 . 180E-02 

0.736E-02 

0 . 983E-03 

0.311E-02 

; 0 . 883E-03 | 

0 .3112-02 

jOHg 

0.736E-02 


Iterate 2 





Row: 4 

0 . 255E-02 

0 . 105E-02 

O'. 431E-03 (0 . 154E-03 {0 . 148E-03 JO . 443E-UO 

0 . 123E-02 

0 . 60 4E-Q2 

Row: 3 

0 . 173E-02 

0 . 109E-02 

0 .143E-03 

0 . 756E-04 



0 .764E-03! 

0.203E-02 

Row: 2 

0 . 173E-02 

0.109E-02 

0.143E-03 

0 . 756E-04 

0.713E-04 i 

0.131E-03I 

4 

0 ,764E-03 ! 

0.203E-02 

Ibb 


0 . 105E-02 

0 .431E-03 

0 . 154E-03 


0 .123E-02- 

0 . 604E-02 

RRORS IN NODALLY f 

lveraged gradient: Iterate 3 

— "■ 


Row: 4 

0 . 185E-02 

10 . 8S9E-03 

1 

0.433E-03 

0 . 142E-03 

0 . 132E-03 

0 . 464E-03 

0 . 954E-03 


Row: 3 

0 . 131E-02 

|o . 859E-03 

0 . 1S6E-03 

0 . 910E-04 

0 . 852E-04 

0 . 144E-03 

0 . 629E-03 1 


n 


|o . 859E-03 

0 . 156E-03 

0 . 910E-04 

0 . 852E-04 

0 . 144E-03 

0 . 629E-03 

EBS3E9 

Row: 1 

0.185E-02 

jo . 859E-03 

0 . 433E-03 

0.142E-03 



0.132E-03 

0 .464E-03 

1 

0.954E-03 

0.556E-02 


Table 9: Element by element errors for iterates in Cantin-Loubignac scheme for the example 
problem. 
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In Table 1 1 we repeat our earlier experiment concerning the various choices of boundary 
condition for the local problem. First, we note that the localized schemes once again are 
superior to the full Cantin-Loubignac scheme. Second, the choice of boundary condition is 
not so critical for the improved recovery scheme determining p. The approximate choice 
of boundary condition is almost as good as the true boundary condition. Furthermore, the 
local behavior of the schemes is satisfactory and identifies the poor approximation pss Vu 
in the neighborhood of the boundaries and corrects accordingly. 

Example: Cracked Panel Problem 

In the examples considered so far, the solution has been smooth. In order to consider 
the effect of singularities, we now study the case of a cracked panel problem: 

Find u : 


— Vu = 0 in ft 
uon Tp 

0 on Tat 

where CI,Td and are shown in Figure 3.70. The true solution to this problem is 


u 

du 

dn 


u(x,y) = r^cos^# 


i.u - /: - -j 


j_i_ . i?- 


so that there :s a □luguiunuj m ncdi tuc u i i xi ^i.c. au t iic i>ip oi Lxie ciacxj. 

The results obtained for the full Cantin-Loubignac scheme are shown in Tables 12-13, 
from which one sees that the scheme has little success in correcting the poor approximation 
to the gradient in the neighborhood of the singularity. 

The behavior of the localized scheme is shown in Table 14. Once again we pay attention 
to the effect of boundary conditions for the local problem with exact boundary conditions 
the localized scheme is superior to the full Cantin-Loubignac scheme. However, the results 
show how sensitive the method is to the boundary condition. The results obtained with the 
approximate choice of fluxes corresponding to % • p proving completely unsatisfactory. 
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Figure 3.70: Domain, mesh and boundary co: 
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for Cracked Panel Problem. 




Iterate 

||Vu - p\\L*(n) 

0 

4.874(-2) 

1 

2.548(-2) 

2 

2.101(-2) 

3 

1.890(-2) 

4 

1.761 (-2) 

5 

1.673(-2) 

6 

1.607(-2) 


1.556(-2) 


1.513(-2) 


Iterate 

Vu - P||l,2(fi) 

0 

.1359 

1 

7.416(-2) 

2 

6.208(-2) 

3 

5.652(-2) 

4 

5.321(-2) 

5 

5.097(-2) 

6 

4.932(-2) 

7 

4.802(-2) 

8 

4.670(-2) 




Table 10: Accuracy of Cantin-Loubignac iterates for the smooth example on finer meshes. 










ERRORS IN DIRECT GRADIENT APPROXIMATION: 


BBS 

rdfllEBfl 

0 . 924E-02 

0.338E-02I 

0 . 452E-03 

0 . 452E-03 

0.338E-02 0 . 924E-02 

0 . 180E-01 


0 . 180E-01 

0 . 924E-02 

0.338E-02 

0 . 452E-03 

0.452E-03 

0.338E-02 

0 . 924E-02 

0 . 180E-01 

Row: 2 

0.180E-01 

0 . 924E-02 


0 . 452E-03 

0.452E-03 

0.338E-02 

0 . 924E-02 

0.180E-01 

Row: 1 

0.180E-01 

0 . 924E-02 

0.338E-02 

0 . 452E-03 

0.4S2E-03 

0 . 338E-02 

0 .924E-02 

0 . 180E-01 

RRORS n 
rror in 

CORRECTS 

H~1 : 9 

D GRADIENT: 

. 1923280705362D-02 

jRow: 4 !0 .216E-04 1 0 .216E-04 

0 .2l6E-04j 0 .216E-04 

0.216E-04 

0 .216E-04 jO .216E-04 

jo .250E-02 

Row: 3 

0 . 911E-05 {0 . 911E-05 

0 . 911E-05 j 0 . 911E-05 

0.911E-05 

0 . 911E-05 

0 . 911E-05 

0 . 151E-02 

tmm w *1 1 hi n 

Row: 2 

lo . 911E-05 |0 . 911E-05 

0 . 911E-05] 0 . 911E-05 

0 . 911E-05 

0 . 911E-05 

0 . 911E-05 10 . 151E-02 j 

(row: 1 1 0 . 216E-04 jo . 216E-04 

0 . 216E-04 j 0 .216E-04 

0.216E-04 

0.216E-04 

0 .216E-04 *0 .250E-02 J 


(a) Results obtained with exact boundary conditions. 


ERRORS IN CORRECTED GRADIENT: 

0.13437981713121 


: Row . 4 0.695E-03 i0.269E-03j0.269E-03 i0.269E-03 j0.269E-03 J 

0.269E-03J 

0.269E-03 |o .274E-02 .; 

j Bow . 3:0 S82E-03 lo .256E— 03f 0 .256E-03 jo .256E-03 {0 . 256E-03 t 

0 .256E-03 

0.256E-03 ;0 . 17 6E-02 t 

Irow: 2:0. 582E-03 { 0 .256E-03 j 0 .256E-03 i 0 .256E-03 jo .256E-03 

0.2S6E-03 

0 .256E-03 jo . 17 6E-02 i 

— i — ~j 

Inow 1 Jo 5 95E-03 1 0 .269E-03! 0 .269E-03 iO .269E-03 ;0 .269E-03 

0.269E-03 

0 .269E-03 *0 .274E-02 ! 

(b) Results obtained using n.p as a boundary condif 

RRORS IN CORRECTED GRADIENT: 
rror in H~l: 0.13437981713121 

;ion. 


IRow: 4 jO . 695E-03 

0 .269E-03 j 0 .269E-03 0.269E-03: 

0.269E-03 jO .269E-03 

1 - 

0.269E-03 

0.274E-02 

Row: 3 i0 . 682E-03 

0 .256E-03| 0 .256E-03 0 .256E-03 i 

0.256E-03 

0.256E-03 



Row: 2 jo . 682E-03 


0.2S6E-03 

0.256E-03 

0.256E-03 

0.176E-02 

i 

Row: l':0.695E-03 

!o .269E-03 j 0 .269E-03? 0 .269E-03 

0.269E-03 

0.269E-03 

1 

0.269E-03 

0.274E-02 

- ■ " — — — — 


(c) Results obtained using (Jj) as a boundary condition. 


Table 11: Impact of different choices of boundary condition for localized problems. 
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Iterate 

||Vu — PWlhq) 

0 

.515909 

1 

.363061 

2 

.336462 

3 

.324060 

4 

.316089 

5 

.310410 

6 

.306194 

7 

.303001 

8 

.300554 

Q 

1 T 1 

OOCAfil 


||Vu- V«|| L 2 (0) = 0.200177 


Table 12: Accuracy of Cantin-Loubignac iterates for the Cracked Panel example. 
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1*1 Cl M (■} t*l H tn M 


ERRORS IN DIRECT GRADIENT APPROXIMATION: 
fr-ror in H^l: 0.20017721413491 



H 





O.490E-O4 

0 . 434E-04 

0 .280E-04 1 

0.117E-04 



0.271E-04 

0 . 617E-04 

0.117E-03 

0 . 150E-03 

0.953E-04 

0 . 833E-04 

0.701E-04 

0.365E-04 


1 

0.236E-04 

0.766E-04 

0.333E-03 

0 . 649E-03 

0 . 414E-03 

0.528E-031 



Row: 

E 

0 . 605E-05 

0.333E-04 

0.329E-03 

0.187E-01 

0 . 168E-01 

0 . 734E-03. 

0 . 211E-03 

0 . 101E-03 


rrors in nodally averaged gradient: Iterate 0 

rror in H~l: 0.1848138 7268082 


Row : 4 


0 .197EH54 


0 . 343E-04 


0 . 6G~2e-04 


Q.839E-Q4 


0.874E-04 


Q.557E-04 


0 . 177E-04 


0 . 957E-05 


Row: 3 


0.330E-05 


0.540E-05 


0.339E-04 


0 . 809E-04 


0.123E-03 


0. 672E-04 


0.230E-04 


0.488E-04 


Row: 2 


0.218E-05 


0 . 570E-05 


0.208E-04 


0 . 637E-03 


0 . 924E-03 


0 . 497E-04 


0 . 954E-04 


0 . 114E-03 


Row: 1 


0.557E-05 


0.383E-04 


0.265E-03 


0 .110E-01 


0.183E-01 


0 . 144E-02 


0.249E-03 


0 . 195E-03 


RRORS IN NODALLY AVERAGED GRADIENT: Iterate 1 








nraaa 

0 . 933E-05 

0.S27E-05 

M 

0 . 162E-04 

0.203E-05 

0.557E-05 


0 . 825E-05 

0 . S15E-05 

0.943E-05 

0 . 2 68E-05 

Row: 2 

0.134E-04 

0 . 577E-05 

0.292E-04 

0 . 542E-03 

0 . 613E-Q3 

0 . 217E-03 

0 . 199E-04 

0 . S98E-05 

Row: li 

0 . 483E-05 

0 . 677E-05 

0 . 907E-04 

0.102E-01 

0.148E-01 

0 . 145E-02 

0 . 796E-04 

0.146E-04 


RRORS IN NODALLY AVERAGED GRADIENT: Iterate 2 


Sum 

IzMe5^HH!0E! 

aaHaaaadsc 

HEBSH 


0.305E-04 

0 . 102E-04 


0.200E-04 

0.387E-05 



0.274E-05 

0 . 380E-05 

0 . 128E-04 

0.109E-04 

0 . 199E-04 

0 . 191E-04 

0 . 561E-05 

Row: 2 

0 . 191E-04 

0 . 147E-04 

0.343E-04 

0 . 484E-03 

0 . 612E-03 

0.363E-03 

0.137E-04 

0 . 903E-05 

Row: 1 

0 .728E-05 

0 . 113E-04 

0.601E-04 

0 . 101E-01 

0 . 137E-01 

0 . 169E-02 

0 . 217E-03 

0 . 981E-05 


IN NODALLY AVERAGED GRADIENT: Iterate 3 


:ror m 
Row: 4 

u.i 

0 . 203E-04 

.00 / OOU4 

0.848E-05 i 

0 . 158E-04 

Q.306E-Q4 

0.113E-04 

0.217E-04 

0.226E-04 

0.911E-05 

Row: 3 

0 . 183E-04 

0.307E-05 

0.580E-05 

0 . 184E-04 

0.253E-04 

0 .342E-04 

0.241E-04 

0 . 693E-05 

Row: 2 

0.207E-04 

0.221E-04 

0.392S-04 

0.450E-03 

0 . 643E-03 

0 . 447E-03 

0.189E-04 

0.1^8E-04 

Row: 1 

0 . 884E-05 

0.184E-04 

0.522E-04 

0 . 102E-01 

0 . 131E-01 

0 . 185E-02 

0 . 325E-03 

0.139E-04 


Table 13: Element by element errors for iterates in Cantin-Loubignac scheme for Cracked 
Panel. 
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Summary and Conclusions 

The basic aim of the work described in this section was to see whether the extremely 
expensive full Cantin-Loubignac scheme could be localized (and therefore be made cheaper) 
without sacrificing too much accuracy. The results presented here show that the answer 
is affirmative for smooth problems. In fact, the localized scheme can be superior (owing 
to the higher order accuracy by which the local problems can be solved). However, the 
determination of boundary conditions for the local problems is critical. The simple schemes 
presented here work well for smooth problem but are unsatisfactory for a non-smooth case. 
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ERRORS IN DIRECT GRADIENT APPROXIMATION: 
Error in H A 1: 0.20017721413491 


Row: 4 

G.200E-04' 

0 .356E-04 

0 . 526E-04~! 

0.569E-04 

0 . 498E-04 

0.434E-04 

0.280E-04 

0.117E-04 

Row: 3 

0.271E-04 

0 . 617E-04 

j 

0.117E-03! 

0.150E-03 

| 0.953E-04 

0 . 833E-04 

0 .701E-04 

0.365E-04 

Row: 2 

0.236E-04 

0.766E-04 

0.333E-03 

0 . 649E-03 

| 0 . 414E-03 

0.528E-03 



Row: 1 

0 . 605E-05 

0.333E-04 

0.329E-03 



0.734E-03 

0.211E-03 

0 . 101E-03 


ERRORS IN CORRECTED GRADIENT: 
Error in H~l: 0.11623575645905 


-Lll 

Row: 4 

0.103E-04 

0.233E-04 

0 . 442E-04» 0 . 661E-04J 0 .291E-04 

0.238E-04 

0 . 118E-04 

0.738E-05 

Row: 3 

0 . 187E-05 

0.525E-05 


0.208E-05 

0.848E-07 

0 .395E-04 

Row: 2 

0.210E-05 

0 . 526E-05 

0.164E-04j0.606E-03 0.152E-03 



0 . 935E-04 

Row: 1 ] 

0 .560E-05 

0 . 380E-04 

0 .258E-03i 0 .107E-0l| 0 . 973E-03 

0 . 149E-03 

0.374E-06 

0.145E-03 


(a) Results obtained using exact fluxes for boundary conditions. 


ERRORS IN CORRECTED GRADIENT: 
Error in H~l: 0.30233586712545 


(Row: 4 10 . 1.03E-Q 4 

0 . 233E-041 0 . 442E-04i 0 . 661E-04 * 0 * 603E-04 r ! 0 . 444E-04 0 . 175E-04 

0.738E-05 

I Row: 3 

0 .187E-05 

0 .S25E-05I 0.334E-0410 . 800E-04 1 0 . 105E-03> 0 . 102E-03 0 . 910E-04 

0 . 395E-04 

Row: 2 

0 .210E-05 

0.526E-05 

0 .l64E-04j 0 . 606E-03 !o . 695E-031 0 . 169E-02 0 . 870E-04 ] 

0.935E-04 

Row: 1 

•0 . 560E-05 1 

0 . 380E-04 

0 .258E-03 ) 0 . 107E-01 |b . 746E-01 \ 0 . 157E-02 j 0 . 141E-03 j 

0 . 145E-03 ' 


(b) Results obtained using n.p as a boundary condition. 



(c) Results obtained using (£) as a boundary condition. 


Table 14: Impact of different choices of boundary condition for localized problems. 
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4 Features of Auditor 


Many of preprocessing and postprocessing techniques outlined in Section 3 above have been 
assembled into a highly specialized software package called A UDITOR. The functionality of 
this package focuses on three distinct aspects of the design and analysis cycle all of which 
play an important part in the success of the modeling effort and the reliability of the results. 
Included in the design of AUDITOR are the following features and capabilities: 


Data Structure 

A UDITOR uses an object based type of data management system. 

Memory Management 

The data structure within AUDITOR uses dynamic memory allocation, thereby placing 
no artificial limitations on the problem size, and minimizing the storage required for smaller 
problems. 

Two and Three-Dimensional Domains 

AUDITOR can provide pre- and postprocessing of both two-dimensional and three- 
dimensional computational domains. 


Problem Classes 


A r\rm r\r\ on 

j » vx - -& VX M. v |/1 \S * j^/i v 


:- and processing for the following problem classes. 


• Steady state heat conduction 

• Linear elasticity 

• Incompressible Viscous Flow 

• Compressible inviscid and viscous flows 


At the present time, the AUDITOR error estimator and solution enhancement module only 
accept laminar flows. 
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Numerical Models 

The interface to AUDITOR has been designed to accept numerical solutions obtained 
from the following classes of solvers. 

• Finite Element 

• Finite Volume 

• Finite Difference 

Note, these models must be presented in a particular format to be read into the module 
which at the present time is based on a hexahedral cell. 

Error Estimation 

AUDITOR features a residual error estimator for each of the problem classes outlined 
above. This error estimator provides for local and global estimates of the error for either a 
single component of the solution or for all of the solution components collectively. 

Solution Enhancement 

A UDITOR contains a postprocessing option for extracting higher order continuous solu- 
tion results for selected derivative components of the solution. 

Mesh Preprocessing 

AUDITOR contains a mesh-preprocessing module for checking the orthogonality and 
smoothness of a computational mesh. If problems are detected with the mesh a summary 
table is provided on the screen and in separate output file. 

Boundary Conditions 

The unique formulation of the residual error estimator within AUDITOR has been de- 
signed to not require that boundary conditions for the problem be supplied. This feature 
basically enforces the assumption the solution to the problem being solved is exact on the 
boundary. 

Interfaces 

AUDITOR provides for two formats to accept input data to be processed. These include 
the COMCO gridfile format and the P3/CFD restart file format. (P3/CFD is a general 
purpose fluid flow analysis package marketed by PDA Engineering.) 

Graphical User Interface 
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A UDITOR was designed to operate in a workstation environment and is essentially mouse 
driven. The graphical interface is motif based and generally can operate in Unix environment 
which has X_windows or GLX libraries. 

Visualization Capabilities 

AUDITOR provides a wide array of graphics visualization capabilities for reviewing the 
solution and/or computational domain. Included in these features are 

• Static and dynamic modes of operation 

• Software and hardware graphics 

• Visualization capabilities 

— isosurfaces 

— cutting surfaces 

— line probes 

— boundary profiles 

— velocity vectors 

A summary of these features is provided in table 1. 
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• Higher order solution enhancement 


Numerical Models 

- 2D or 3D 

- linear or quadratic solutions for hexahe- 
dral meshes* 

— Finite Volume 

- Finite Difference 

- Finite Element 

* accepts HP meshes from P3/CFD 
or PH LEX 

Error Estimation Capabilities 

— Residual based error estimator for 

* Heat Conduction 

* Linear Elasticity 

* Incompressible Flow (Laminar) 

* Compressible Flow (Laminar) 

- Local cell-wise estimate of the error 

- Global estimate of the solution error 

Mesh Preconditioning 

- Identifies poorly conditioned regions of 
the mesh 

— Identifies tangled meshes 

— Provides a summary file of the precon- 
ditioning analysis 


- provides higher order results for various 
derived quantities 

• Neutral File 

- Reads COMCO predefined format (see 
section 3) 

- Reads P3/CFD restart file 

- Reads PHLEX restart file 

• Graphical User Interface 

- Completely mouse driven 

- Operates in XT mode or GLX on an SGI 
workstations 

• Graphics/Visualization 

- Isosurfaces 

- Velocity Vectors 

- Postscript Hardcopy 

- Slices 

- XY Plotting 

- Visualization Variables include: 

* solution components 

* predefined derivatives and combina- 
tions 

* velocity magnitude 

* local error estimate 

* elements with bad jacobians 


Table 15 : AUDITOR Features and Capabilities 




Figure 5.1: Isotropic insulated bar geometry and boundary conditions for a steady state heat 
conduction analysis. 

5 Test Problems 

5.1 Introduction 

This section of the final report presents selected results for a number of example problems 
which have been analyzed with the A UDITOR package. The scope of these problems range 
from two-dimensional steady state heat conduction to three-dimensional compressible vis- 
cous flows. The solutions to the heat conduction and linear elasticity problems are all finite 
element based solutions which have been obtained by the PHLEX code. This is a finite ele- 
ment code developed by the Computational Mechanics Company which includes hp-adaptive 
solution capabilities and also disposes a compatible neutral file for supplying solution data to 
AUDITOR. The solutions for the incompressible viscous flows have come from two sources, 
P3/CFD, a general purpose finite element CFD code marketed by PDA Engineering, and 
from the TEACH code developed at the Imperial College of Science and Technology in Lon- 
don. Finally, the solutions for the compressible flows have been taken from P3/CFD, and 
from a finite volume code, JCODE, which is based on Jameson’s finite volume Runge-Kutta 
time stepping scheme. 

On a final note before beginning a discussion of the test problems; the graphics and 
visualization capabilities in the AUDITOR package are based on a three-dimensional model 
and thus even for two-dimensional solutions the graphics will appear as extruded three- 
dimensional plots with a unit thickness. 

5.2 Steady State Heat Conduction 
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Figure 5.2: Isotropic insulated bar, finite element model consisting of 8 linear elements. 

5.2.1 Isotropic Insulated Bar (PHLEX - finite element solution) 

The first sample problem presented herein is a steady heat conduction analysis of two- 
dimensional rectangular bar. The geometry and boundary conditions used to solve this 
problem with the PHLEX code are shown in Figure 5.1. Summarizing the boundary condi- 
tion data shown in this figure, there is a prescribed temperature of 100° on the left face, a 
prescribed temperature of 1000° on the right face, and insulated boundaries (zero heat flux) 
on the upper and lower faces. The thermal conductivity is assumed to be independent of the 
coordinate location and is arbitrarily selected to be 1.0. For this problem, an exact solution 
is known and is given by the following expression 

u(x, y, z) = 100 + 900 * xjl 

Here x is the axial coordinate length measured from the left end of the bar and / is the 
length of the bar. 

Numerical Model 

The geometry, shown in Figure 5.1, has initially been discretized using 8 linear finite 
elements with 4 elements in the axial direction and two elements through the thickness, see 
Figure 5.2. 

AUDITOR Analysis 

This problem was solved using the PHLEX code which generates a neutral file format 
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which can be directly read into AUDITOR. As this mesh is quite simple, there is basically 
no need to preprocess the mesh so this phase of the analysis will be omitted. The features 
of interest for this linear steady state problem are the temperature distribution and the 
global and local error distribution. To ensure that the numerical has been correctly read 
into the database, an isosurface plot of the temperature distribution is shown in Figure 
5.3. (Again, note that the graphics package is a three-dimensional package which shows 
the computational domain with the geometry and solutions expand with a unit value in the 
third dimension.) From this figure we observe that the temperature variation in the mesh 
has exactly captured the linear variation of the temperature from 100 to 1000. Activating 
the residual error estimator, the estimated error in the solution is zero verifying the residual 
error estimator in operating correctly. 
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Figure 5.3: Temperature isosurfaces for the insulated bar. 
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5.2.2 Isotropic Bar With Prescribe Source (PHLEX - finite element solution) 

A second example based on a simple rectangular geometry, shown in Figure 5.4, has been 
used in the solution of a three-dimensional heat conduction problem with a slightly different 
set of boundary conditions and a heat source term. Summarizing the boundary condition 
data shown in this figure, all of the boundaries of the three-dimensional bar have a pre- 
scribed temperature of zero and the heat source term has been selected to have the following 
quadratic volumetric distribution 

s(x,y,z) = -310 *((-2)* (1 - y) * y(3 - z) * z 

+(2 — x) * x * (—2) * (3 — z) * z 

+(2 - x) * x * (1 - y) * y * (-2)) 

—170 * ( (4 — 6 * x) * (1 — y) * y * (3 — z) * z 
+(2 — x) * x * *2 * (— 2) * (3 — z) * z 

+ (2 — x) * x * *2 * (1 — y) * y * (—2)) 

—290 * ( (—2) * (1 — y) * y * *2 * (3 — z) * z 

4- (2 — x) * x * (2 — 6 * y) * (3 — z) * z 

+(2 — a:) * x * (1 — y) * y * *2 * (—2)) 

—130 * ( (—2) * (1 — y) * y * (3 — z) * z * *2 
+(2 — x) * x * (—2) * (3 — z) * z * *2 

+(2 — x)*x*(l — — Q* z )) 

—230 * ( (—2) * (1 — y) * y * *2 * (3 — z) * z * *2 
+(2 — x) * x * (2 — 6* y)*) * (3 — z) * z * *2 
+(2 — x) * x * (1 — y) * y * *2 * (6 — 6 * z)) 

—110 * ( (4 — 6 * x) * (1 — y) * y * (3 — z) * z * *2 
-f(2 — x) * x * *2 * (—2) * (3 — z) * z * *2 

+(2 — x) * x * *2 * (1 — y) * y * (6 — 6 * z)) 

—190 * ( (4 — 6 * x) * (1 — y) * y * *2 * (3 — z) * z 
+(2 — x)*x**2* (2 — 6*?/)* (3 — z) * z 

+(2 — a;) * x * *2 * (1 — ?/)*?/* *2 * (—2)) 

—100 * ( (4 — 6 * x) * (1 — y) * y * *2 * (3 — z) * z * *2 
+(2 — x) * x * *2 * (2 — 6 * y) * (3 — z) * z * *2 
+(2 — x) * x * *2 * (1 — y) * y * *2 * (6 — 6 * z)) 

The selection of this source distribution allows an exact solution to be determined which 
is given by the following expression: 
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Figure 5.4: Isotropic bar with prescribed volumetric source, geometry and boundary condi- 
tions. 
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u(x, y, Z ) = -310 * ( (2 - x) * x * (1 -y)*y * (3 -z)*z) 

_170 * ( (2 - x) * x * (1 - y) * y * (3 - z)*z)*x 

-290 * ( (2 - x) * x * (1 -y)*y * (3 -z)*z)*y 

-130 * ( (2 - x) * x * (1 -y)*y * (3 - z) * z) * z 

-230 * ( (2 - x) * x * (1 -y)*y * (3 -z)*z)*yz 

— 110 * ( (2 — x) * x * (1 -y)*y *(3 -z)*z)*xz 

-190 * ( (2 - x) * x * (1 - y) * y * (3 -z)*z)*xy 

-100 *((2-x)*x * (1 -y)*y * (3 - z) * z) * xyz 

The thermal conductivity is again assumed to be independent of the coordinate locations 
and we have selected it to be 1.0. Note, to run this problem usmg AUDITOR, an environment 
variable which triggers this body force needs to be set as follows: 

setenv POISSON-TEST ON 


Numerical Model 

The geometry, shown in Figure 5.4, has initially been discretized using 16 finite elements 
with 4 elements in the axial (x) direction and 4 elements in the cross section. The PHLEX 
analysis for this problem was carried out in three stages which consisted of calculating the 
solution on a linear finite element mesh, an isotropic quadratic mesh, and an isotropic cubic 
mesh. The goal with the AUDITOR analysis is to observe the change in the residual error 
as the spectral order of the mesh is increased. Note, that since a very specific volumetric 
source term was selected, the exact solution is a cubic polynomial which implies that the 
AUDITOR code should predict a zero error for the third order PHLEX analysis. 

A UDITOR Analysis 

The meshes used in all three of the solution sequences are again quite simple and no 
preprocessing of the mesh is required. The features of interest for this steady state analysis 
include the temperature distribution within the domain, and a monitoring of the global error 
as the spectral order of the mesh is increased. To ensure that the numerical solutions have 
been correctly read into the database, a plot of the isosurfaces of the temperature of each 
of the models is shown in Figure 5.5 - 5.7. Calculating the residual errors for each of the 
meshes, the error in the linear solution is 261%, the error in the quadratic mesh is 8.3 %, 
and the error in the cubic solution is zero. 

5.2.3 Connecting Rod (PHLEX - finite element solution) 

The final steady heat conduction analysis presented here is of a three-dimensional connecting 
rod. The geometry and boundary conditions used in the PHLEX simulation are shown in 
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Figure 5.8. Summarizing the boundary conditions, there is a prescribed temperature of 
500° on the piston pin contact surface, left most circular arc in Figure 5.8, a prescribed 
temperature of 250° on the crank shaft contact surface, right most quarter circular arc, and 
zero normal heat flux on the remaining boundaries. The volumetric source term is set to 
zero and the thermal conductivity is specified to be a constant of 1.0. 

Numerical Model 

The geometry, shown in Figure 5.8, has initially been discretized using 194 linear finite 
elements which have been enriched to second order for the solution phase, see Figure 5.9. 

AUDITOR Analysis 

The first phase of the AUDITOR analysis for this test problem accesses the AUDITOR 
pre-processor to gather diagnostics concerning the mesh quality. Figure 5.10 shows a copy 
of the diagnostics which appear on the screen as a result of this preprocessing phase. Of 
note here are the following: there are no elements with bad surface or volume jacobians, the 
critical angle selected for orthonogality testing was 30° and 8 cells in the mesh were found to 
have angles which are outside the critical boundary of 30° and 150°, the largest and smallest 
included angles in the mesh are 177.19 and 15.87 respectively, and the largest volumetric 
and surface Jacobian ration was found to be 0.069503. This final value, the jacobian ration, 
in general, gives an indication that at least one of the elements used in the computations is 
severely distorted locally and could possibly cause numerical problems. 

In the second phase of the analysis, the global and local error distributions were cal- 
culated. Based on the quadratic finite element approximation to the solution, the residual 
error estimator indicates that there is a 1.56% global error in the solution with a maximum 
local error of 0.52%. Figure 5.11 shows the corresponding distribution of the residual error 
projected onto the surface of the connecting rod. Interpreting the color coding in this figure, 
it appears that the highest errors in the solution occur near the intersections of the connect- 
ing rod and the Piston contact surfaces. A zoom of this region is shown in Figure 5.12. For 
reference, we have also included a copy of the temperature distribution, Figure 5.13 which 
shows the maximum temperature of 500° on the piston contact surface, and a temperature 
of 250° on the crank shaft contact surface. 
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Figure 5.5: Temperature isosurface for an isotropic bar with a prescribed source for a linear 
finite element solution 
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Figure 5.6: Temperature isosurface for an isotropic bar with a prescribed source for a 
quadratic finite element solution 
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Figure 5.7: Temperature isosurface for an isotropic bar with a prescribed source for cubic 
finite element solution 
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Figure 5.8: Three-dimensional connecting rod, geometry and boundary conditions 



Figure 5.9: Finite element for a connecting rod mesh consisting of 194 elements and 2417 
degrees of freedom 


170 


«< MESH CONDITIONING SUMMARY »> 


* SUMMARY OF VOLUMETRIC JACOBIAN TESTS 

* Number Cells with Bad Volume Jacobians 

* Maximum Volumetric Jacobian ratio 

* - Detected in Cell 


* SUMMARY OF SURFACE JACOBIAN TESTS 

* Number Cells with Bad Surface Jacobians 

* Maximum Surface Jacobian ratio 

* - Detected in Cell number 


0 

0.069503 
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0 

0.069503 
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* SUMMARY OF ORTHOGONALITY TESTS 

* Critical Minimum Angle Specified : 30.000000 

* Critical Maximum Angle Specified : 150.000000 

* Number of the Cells failing Critical Angle test : 8 

* Largest include angle detected * 177.196407 

* Smallest include angle detected : 15.877637 


* ! ! BAD JACOBIAN OR ANGLE FOUND!! See f lie : . /MeshCond_Results 

* 

* 

* NOTE: CELL NUMBERS REPORTED HERE FOR RESTART FILES MAY EXCEED 

* THE NUMBER OF ACTIVE CELLS REPORTED IN THE Info FORM 


«< END OF SUMMARY »> 


Figure 5.10: Preprocessor analysis 


of the computational mesh for a connecting rod 
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Figure 5.11: Residual error estimation distribution for a steady state solution of a 

three-dimensional connecting rod 
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Figure 5.12: Zoom of the error estimation distribution in the piston contact region 
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5.3 Linear Elasticity 

5.3.1 Isotropic Bar (PHLEX - finite element solution) 

The first linear elasticity problem presented herein is a three-dimensional rectangular bar 
in simple tension. The geometry and boundary conditions used to solve this problem with 
the PHLEX code are shown in Figure 5.14. Summarizing the boundary conditions, there 
are no penetration boundaries along the left, bottom, and back surfaces, a uniform surface 
traction 1000 lbs. in the horizontal x direction on the right face and zero surface tractions on 
the remaining surfaces. The body force terms are zero and the following material constants 
apply; 


A = Lamae’s constant = 5,617,280 psi 
G= Shear modulus = 2,407,407 psi 

The exact solution is easily calculated to be a linear function of the coordinate locations 
with the following components; 

u(x,y,z) = 9.230770e -4 * x/l (in.) 

u(x,y,z) =- 5384516e _s * y/l (in.) 
w(x,y,z) = — 5.384516e~ 5 * z /l (in.) 

Here, x is the axial coordinate location along the bar, y and z are the transverse coordinate 
directions, and l is the length of the length of the bar. 

Numerical Model 

The geometry for this example, shown in Figure 5.14, has initially been discretized using 
16 linear finite elements with 4 elements in the axial direction and four elements in the cross 
section, see Figure 5.15. 

AUDITOR Analysis 

Due to the simplicity of the computational mesh, and solution, the only phases in the 
A UDITOR analysis that are of interest are the error estimation phase and possibly the solu- 
tion enhancement phase. Before activating either of these options, a plot of the isosurfaces 
for the u-displacement component have been extracted as shown in Figure 5.16. This figure 
shows that the axial displacements vary linearly from 0.0 inches on the left end to 9.23077 e -4 
(in.) at the right end, which corresponds to the exact solution given above. Activating the 
residual error estimator we find that the estimated error in the solution is zero verifying that 
residual error estimator is operating correctly. For this particular case, the stress distribution 
is constant, a xx = 1000 psi (all other <r,j = 0.0) and thus initiation of the postprocessing 
solution enhancement shows no improvement in any of the stress components. 
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Figure 5.13: Surface temperature distribution for an isotropic connecting rod 
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Figure 5.15; Finite element mesh for a linear elastic bar under simple tension 



Figure 5.16: U-displacement isosurfaces for a rectangular bar in simple tension 
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5.3.2 Wheel Rim (PHLEX - finite element solution) 

The final linear elasticity problem presented herein is for a three-dimensional analysis of 
wheel rim on a Pontiac Fiero. The geometry and boundary conditions used in the PHLEX 
simulation are shown in Figures 5.17 and 5.18. Summarizing the boundary conditions in 
these figures, there is a prescribed lateral bead pressure of 313 psi due to tire pressurization, 
a vertical load of 123.4 psi on the rim due to weight of the automobile, fixed displacements 
at the lug-nut locations, and free surfaces elsewhere. The body force terms are prescribed 
as zero and we have selected the following material constants: 

A = Lamae’s constant = 5,617,280 psi 

G= Shear modulus = 2,407,407 psi 


Numerical Model 

The geometry, shown in Figures 5.17 and 5.18, has initially been discretized using 570 
linear finite elements with 1300 degrees of freedom, see Figure 5.19. Note, that the model 
shown in Figure 5.19 contains symmetry only about the z = 0 and y — 0 planes. 

AUDITOR Analysis 

The first phase of the A UDITOR analysis for this test problem accesses the preprocessor 
to gather diagnostics concerning the mesh quality. Figure 5.20 shows a copy of the diag- 
nostics which appear on the screen as a results of this preprocessing phase. Of note here 
f^U/vurirwr- tKere »rp no piemen t.s with bad surface or volume iacobians, the largest 
and smallest included angles in the mesh are 136.0° and 25.5° respectively, and the larges 
jacobian ratio was found to be 0.228 which indicates only moderate distortions of mesh and 
which should not effect the computations. 

In the second phase of the analysis we have calculated the global and local error distribu- 
tions. For this linear finite element solution the residual error estimator indicates that there 
is an 532% global error in the solution with a maximum local error of 75.5%. Figure 5.21 
shows the corresponding distribution of the residual error projected onto the wheel surface. 
Based on the color coding of this figure, the highest errors in the solution are occurring 
almost directly below the wheel axle at the intersection of the cylinder drum and the tie bar. 
Note, that the error distribution is symmetric about the z = 0 plane but exhibits no other 
symmetry due to the loading pattern and orientation of the tie bars. 

To complete the analysis of the wheel rim, we have accessed the postprocessing options 
to extract higher order results for the Von Mises stress. As a reference point we have 
plotted isosurfaces on the boundary of the current distribution of the Von Mises stress as 
shown in Figure 5.22. Note, from this figure an the diagnostics in the isosurface plotting 
window, the stresses are discontinuous across element interfaces, and the maximum stress 
value is’ 2520 psi. Figure 5.23 shows the corresponding Von Mises stress distribution after 


177 


the postprocessor has been called. Based on the enhancement procedure, the Von Mises 
stresses are now continuous and the maximum stress value has been slightly relieved. 
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u = 0. 
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Figure 5.18: Geometry and loading configuration an automobile wheel subjected to static 
tire pressure and vertical bead loading. Front view showing radial the distribution of the 
vertical bead load. Tire pressure loads are not shown as they are uniformly distributed 
around the rim of the wheel. 
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Figure 5.19: Linear finite element used to model three-dimensional wheel rim on a Pontiac 
Fiero 
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*********************************** 




«< MESH CONDITIONING SUMMARY »> 


* SUMMARY OF VOLUMETRIC JACOBIAN TESTS 

* Number Cells with Bad Volume Jacobians 

* Maximum Volumetric Jacobian ratio 

* - Detected in Cell 


0.228694 

71 


SUMMARY OF SURFACE JACOBIAN TESTS 

Number Cells with Bad Surface Jacobians 
Maximum Surface Jacobian ratio 
- Detected in Cell number 


0.228694 

71 


* SUMMARY OF ORTHOGONALITY TESTS 

* Critical Minimum Angle Specified 

* Critical Maximum Angle Specified 

* Number of the Cells failing Critical Angle test 

* Largest include angle detected 

* Smallest include angle detected 


15.000000 

165.000000 
0 

136.047060 

25.523589 


NOTE: CELL NUMBERS REPORTED HERE FOR RESTART FILES MAY EXCEED 

THE NUMBER OF ACTIVE CELLS REPORTED IN THE Info FORM 


«< END OF SUMMARY »> 


Figure 5.20: Preprocessing diagnostics for a three-dimensional wheel rim. 
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Figure 5.21: Residual error distibution projected on the surface for a linear elasticity analysis 
of a three-dimensional wheel rim 
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Figure 5.22: Discontinuous Von Mises stress isosurfaces for the three-dimensional wheel rim 
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5.4 Incompressible Viscous Flow 

5.4.1 Backstep Channel (TEACH - finite difference solution) 

The first incompressible viscous flow problem presented herein is a two-dimensional backstep 
channel problem. The geometry and boundary conditions used in the TEACH code simula- 
tion are shown in Figure 5.24. Summarizing the boundary conditions in this figure, there is 
a prescribed quadratic velocity profile along the upper half of the left hand face given by 

u(x, y, z) = 1.5 * [1.0 - (( y - 1.5)/0.5) * *2], 
an outflow boundary along the right face with 

p = 0.0, 

du dv 
dx dx 

and solid wall boundaries elsewhere. The Reynolds number used in the solution process 
was set to 100 based on the height of the channel. 

Numerical Model 

The geometry, shown in Figure 5.24, has initially discretized into a 35 x 31 finite difference 
mesh which has been clustered in horizontal direction near the inflow region and around the 
separation point. The distribution of nodal points in the vertical direction is uniform with 
the first point occurring at a distance of 0.0667 from the solid wall, see figure 5.25. Note 
that this geometry the entrance length up to the backstep has been removed and a fully 
developed quadratic profile has been imposed to simulate the upstream effects. 

A UDITOR Analysis 

The two-dimensional grid used during the solution phase is again a very simple rectan- 
gular grid which has some moderate clustering, and thus there is no need to preprocess the 
mesh. The features of interest include, the recirculation zone in the lower left portion of the 
mesh, the location of the separation point, and the distribution of the error in the numerical 
solution. To ensure that the numerical solution has been correctly read and placed in the 
database, a plot of the isosurfaces for the u - velocity and pressure are shown in Figures 5.26 
and 5.27. From a plot of the u - velocity for a value of -0.004, and with the aid of the graphical 
picker, the separation point for this model is estimated to be 2.875 which corresponds very 
well with numerically accepted value of 3.0. 

Activating the residual error estimator, the estimated global error in the solution is 
117% with a maximum local error of 24.6%. Figures 5.28 and 5.29 show the corresponding 
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Figure 5.23: Enhanced Von Mises stress isosurfaces for a three-dimensional wheel rim 
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Figure 5.24: Backstep geometry used with the TEACH flow solver, Re=100. 

distribution of the residual error projected on the surface and a zoom of the error near the 
inflow region. Based on the coloring in Figure 5.28 the largest cell-wise error are all occurring 
at or near the inflow boundary. The zoom of the errors in this region, rescaled in Figure 
5.29 to a maximum error of 0.10, shows additional relatively large errors occurring directly 
downstream of the inflow region and bounding the recirculation zone. 

The final phase in the analysis, consists of accessing the postprocessing option to extract 
higher order results for the vorticity magnitude. For reference, the current distribution 
of the vorticity magnitude have been calculated as shown in Figure 5.30. Note, that the 
vorticities are highly discontinuous and the isosurface window indicates that the maximum 
vorticity is 6.75. Figure 5.31 shows the corresponding vorticity distribution after the solution 
enhancement has been activated. Note, from the isosurface window the vorticity magnitude 
now has a maximum vorticity magnitude of 7.07 and small undershoot as well. 
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Figure 5.25: Computational mesh used in the TEACH code simulation and AUDITOR 
analysis of a two-dimensional backstep, Re = 100 
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- 0.2 


Figure 5.26: U- velocity isosurfaces for the two-dimensional backstep, Re = 100 
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Figure 5.27: Pressure isosurfaces for the two-dimensional backstep, Re = 100 
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Figure 5.28: Residual error distribution for the two-dimensional backstep, Re = 100 
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Figure 5.29: Zoom of the residual error near the inflow for the two-dimensional backstep, Re 

= 100 
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Figure 5.30: Vorticity isosurfaces for the two-dimensional backstep channel before solution 
enhancement 
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Figure 5.31: A sketch of the flow domain for the flow through a square-duct with a bend, 
Re = 790 
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5.4.2 90° Bend (P3/CFD - finite element solution) 

The final incompressible viscous flow problem presented herein is a three-dimensional flow 
through a rectangular duct with a 90° bend. The geometry and boundary conditions used 
in the P3/CFD analysis are shown in Figure 5.32. Summarizing the boundary conditions 
shown in this figure, there is a prescribed uniform velocity profile along the vertical face on 
the far right, 

u = 
v = 
w = 

an outflow boundary on the horizontal face at 

p = 0.0 (232) 

and solid wall boundaries elsewhere. The Reynolds number used in the solution process was 
set to 790 based on the height of the duct. 


- 1.0 

0.0 (231) 

0.0 

the top with 


Numerical Model 

The geometry for this example, shown in Figure 5.32, has been discretized into 28 
quadratic finite elements with 375 degrees of freedom. In this model, there are 7 elements 
along the length of the duct which have been clustered in the bend region and four elements 
in the cross section, see Figure 5.33. Note, this is a very coarse model for this problem, and 
that the goal here is not to obtain a highly accurate solution from the solver but estimate 
the error in the numerical solution. 

A UDITOR Analysis 

The three-dimensional grid used during the solution phase is again a very simple rect- 
angular grid, which has some moderate clustering, and thus there is basically no need to 
preprocess the mesh. The features of interest include the development of vortex tubes in 
the bend region of the duct, a possible recirculation zone at the end of the bend, and the 
distribution of the error in the numerical solution. To ensure that the solution has been 
correctly read and placed in the database, a cross section plot of the velocity magnitude and 
isosurfaces of the pressure are shown in Figures 5.34 and 5.35. For reference, we have also 
plotted the z-velocity components, see Figure 5.36, which shows the out of plane recirculation 
in the bend region. 

Activating the residual error estimator, the estimated global error in the solution is 
778% with a maximum local error of 256%. Figure 5.37 shows a zoom of the corresponding 
distribution of the residual error projected on the surface of the duct and Figure 5.38 shows 
the error projected onto a planar slice through the center of the duct. Based on the gray 


195 




scale coloring in this figure, the maximum cell-wise errors are occurring at the inflow and 
the outflow sections of the domain. The error at the inflow is most likely a result of the 
technique which is being used to enforce uniform inflow velocity profile. The relatively large 
error at the outflow section is at least in part due to the relatively large size of the quadratic 
elements used near the outflow boundary. A finer mesh in either of these regions would most 
probably reduce the errors there significantly. Finally, note that due to the coarseness of 
the finite element model there is not separation of the flow at the exit point of the bend, 
however, there still appears to be relatively large local errors in this region. 






I 

Figure 5.34: Planar slice showing the velocity magnitude at the centerline of the 90° bend 
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Figure 5.35: Pressure isosurfaces in the bend region of the rectangular duct 
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Figure 5.36: Z-velocity isosurfaces in the bend region of the rectangular duct showing 
three-dimensional recirculation 
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Figure 5.37: Zoom of the residual error distribution in the bend region of the duct projected 
on the surface 
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Figure 5.38: Planar slice through the center of the rectangular duct colored by the error 
distribution 
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The final phase in the analysis of this problem consists of accessing the postprocessing 
options to extract higher order results for the vorticity magnitude. Figure 5.39 shows the 
corresponding vorticity magnitude distribution after the postprocessing phase. Note, the 
scale in this figure shows the enhancement phase has caused an undershoot in the vorticity 
magnitude. 


5.5 Compressible Flow 

5.5.1 Airfoil (JCODE - finite volume solution) 

The first compressible flow problem presented herein is a two-dimensional inviscid flow past 
a NACA 0012 airfoil. The geometry and boundary conditions used in the JCODE simulation 
are shown in Figure 5.40. Summarizing the boundary conditions shown in this figure; far 
field data are prescribed along the entire outer parameter of 0-grid, and solid wall conditions 
(no penetration) are specified on airfoil. The far field conditions are 



Figure 5.39: Isosurface of the vorticity magnitude after postprocessing 
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Mach = 0.75 
alpha = 2.0° 
gamma = 1.4 

All other quantities used in the computations are non-dimensionalized by the density and 
velocity at infinity. 

Numerical Model 

The geometry shown in Figure 5.40, has been discretized into 1200 cells which are clus- 
tered near the airfoil in the radial direction and near the leading and trailing edge in the 
primary flow direction. A zoom of the mesh near the airfoil is shown in Figure 5.41. 

AUDITOR Analysis 

The first phase of the AUDITOR analysis for this test problem accesses preprocessor to 
gather diagnostics concerning the mesh quality. Figure 5.42 shows a copy of the diagnostics 
which appear on the screen as of result of this preprocessing phase. Of note here are the 
following: there are no elements with bad surface or volume jacobians, the largest and 
smallest included angles in the mesh are 171.9 andl4.4 respectively, and the largest jacobian 
ratio was found the be 0.175. 

The flow features of interest for this case are the Mach number distribution and reat- 
tachment point of the shock on the upper surface, the pressure distribution along the surface 
of the airfoil, and the distribution of the error in the numerical solution. To ensure that 
the numerical solution has been correctly read and placed in the database, a plot of the 
isosurfaces for the Mach number and pressures are shown in Figures 5.43 and 5.44. Based 
on the data contained in Figure 5.43, and the isosurface editor form, the maximum Mach 
number predicted by the JCODE is 1.29 which compares well with other results reported in 
the literature. For this case, we have also plotted the pressure distribution along the surface 
of the wing, see Figure 5.45. 

Activating the residual error estimator, the estimated global error in the solution is 
0.29% with a maximum local error of 0.106%. Figure 5.46 shows a zoom of the residual error 
distribution near the airfoil. Based on the color scale in this figure, the maximum cell-wise 
errors are occurring near the leading edge of the airfoil on both the section and pressure 
surfaces and at the shock reattachment point near the center of the airfoil. 
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Figure 5.41: Computational mesh (near the airfoil) used by the JCODE to model inviscid 
M = 0.75 over a NACA 0012 geometry 
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«< MESH CONDITIONING SUMMARY »> 


* SUMMARY OF VOLUMETRIC JACOBIAN TESTS 

* Number Cells with Bad Volume Jacobians 

* Maximum Volumetric Jacobian ratio 

* - Detected in Cell 


0 

0.174991 

1181 


* SUMMARY OF SURFACE JACOBIAN TESTS 

* Number Cells with Bad Surface Jacobians 

* Maximum Surface Jacobian ratio 

* - Detected in Cell number 


0 

0.174991 

1181 


* SUMMARY OF ORTHOGONALITY TESTS 

* Critical Minimum Angle Specified 

* Critical Maximum Angle Specified 

* Number of the Cells failing Critical Angle test 

* Largest include angle detected 

* Smallest include angle detected 


* !!BAD JACOBIAN OR ANGLE FOUND!! See f ile : . /MeshCond Results 


15.000000 

165.000000 

4 

171.928621 

14.447908 


NOTE: CELL NUMBERS REPORTED HERE FOR RESTART FILES MAY EXCEED 

THE NUMBER OF ACTIVE CELLS REPORTED IN THE Info FORM 


«< END OF SUMMARY »> 


Figure 5.42: Preprocessor statistics for the NACA 0012 finite difference grid 
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Figure 5.43: Mach number isosurfaces for NACA 0012 airfoil 
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Figure 5.44: Pressure isosurfaces for NACA 0012 airfoil 
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Figure 5.45: Pressure distribution along the wi 
JCODE. 
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5.5.2 Carter Flat Plate (P3/CFD - finite element solution) 

The final compressible flow problem presented herein is a two-dimensional viscous flow past 
a flat plate. The geometry and boundary conditions used with the P3/CFD flow solver are 
shown in Figure 5.47. Summarizing the boundary conditions shown in this figure there is 
a supersonic inflow along the left face, openflow along the top face, outflow along the right 
face, a symmetry entrance length along part of the bottom face, and a solid wall plate along 
the remainder of the bottom of the domain. 

On the inflow and openflow portions of the boundary the initial non-dimensional flow 
conditions are specified to be; 


Mach 

= 

3.0 

a 

= 

0.0 degrees 

P 

= 

1.0 

u 

= 

1.0 

V 

= 

0.0 

w 

= 

0.0 

Re 

= 

1000.0 

Pr 

= 

0.72 

oo 

— 

390 R 


The conditions on the solid wall plate were taken to be isothermal with a specified 
temperature of 


Twall = 1092 R 

All quantities used in the computations are non-dimensionalized by the density and 
velocity at infinity. 

Numerical Model 

The geometry for this example, shown in Figure 5.47, has initially been discretized into 
77 linear finite elements with 7 elements in the cross stream direction and 1 1 elements in the 
streamwise direction. This mesh was manually edited with a uniform refinement along the 
solid wall boundary, and then two global isotropic refinements were performed. The final 
mesh shown in Figure 5.48 contains 1712 elements and 1829 degrees of freedom. 
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Figure 5 . 46 : Zoom of the residual error estimation distribution near the surface of the wing 
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AUDITOR Analysis 

The flow features of interest for this case include: the shock that develops at the tip of 
the plate and exits the right outflow boundary, the thickness of the boundary later region, 
and the distribution of the error in the numerical solution. To ensure that the numerical 
solution has been correctly read and placed in the database, a plot of the isosurfaces for the 
Mach number and the pressure along the plate are shown in Figures 5.49 and 5.50. These 
figures clearly show the shock emanating from the leading edge of the plate and exiting out 
the outflow boundary. 

Activating the residual error estimator, the estimated global error in the solution is 21.7% 
with a maximum local error of 14.2%. Figure 5.51 shows the corresponding distribution of 
the residual error with the maximum local value cutoff at 5%. Based on the gray scale 
coloring in this figure, the maximum element errors are occurring at the leading edge of 
the plate and in the shock region. One can also observe some relatively large errors in the 
boundary layer region just down stream of thes singularity. 
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Open flow 



Figure 5.47: Carter’s flat plate problem. Geometry and boundary conditions. 
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Figure 5.48: Computational mesh for the Carter flat plate problem, 1712 elements, 1829 
degrees of freedom 
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Figure 5.49: Mach isosurfaces for the Carter flat plate, Re = 1000 
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Figure 5.51: Residual error estimation distribution for the Carter flat plate with the maxi- 
mum cell error limited to 5% 
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6 Future Directions 


The Phase I feasibility study and the Phase II research and development effort have resulted 
in the development of a two- and three-dimensional pre- and postprocessing package called 
Auditor. The scope of this design includes a collection of new and unique features which 
we believe have advanced the state-of-the-art in computational mechanics. As with all 
research and development projects, several areas of research and development exist for which 
additional effort is needed to improve the robustness and utility of the current version of the 
Auditor code. 

In particular, these areas include: 

New Data Structure. 

As indicated in Section 2, we have completed the design phase and much of the coding 
effort on a new object based data structure which is considerably more flexible than current 
object based hexahedral data structure. A reasonably focused development effort in this area 
could significantly expand the class of models (hexahedral, tetrahedral, prism, etc.) which 
Auditor accepts. This would obviously allow a much larger group of solutions to be read in 
and analyzed with Auditor. 

Neutral File Reader. 

As was also indicated in Section 2, we have completed the design phase and much of the 
coding effort on a neutral grid file reader. When completed, this reader will greatly simplify 
the conversion process required to read in finite volume and finite difference solutions as well 
as provide a means for reading in shell elements for elasticity problems. 

Parallel Computations. 

The current version of Auditor is designed to run in a workstation environment with some 
preliminary grouping and coloring of the mesh provided to aid in parallel computations. We 
believe that a moderate level of effort in this area could provide a workstation based version 
of Auditor which could operate with an efficiency of between 3 to 3.75 on a four processor 
workstation. 
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Batch Version. 

Auditor was designed to operate interactively in a workstation environment. For large 
problems, the local computations of the element error and the solution enhancement can be 
computationally intensive. A batch version of Auditor which operates in a CRAY environ- 
ment could provide significant increase in turn around time. Such an effort would require 
not only designing and implementing a batch interface for the code but also would require a 
reasonable effort toward vectorizing the package. 

Error Estimation Scaling. 

Auditor currently uses an estimate of the error based on residual calculation which is 
scaled by the energy norm of the solution for linear problems and an L2 norm of the solution 
for fluid flow problems. To provide the user a better insight on how various estimates of the 
error should be interpreted, further studies into error scaling are needed. 

Directional Error Estimation. 

Auditor presently provides a local estimate of the error on a cell-wise basis. The next 
obvious extension to this isotropic error estimate is to also develop an error estimate which 
has directional characteristics. Such a research effort could result in an error estimate that 
would be quite valuable for doing anisotropic mesh modifications. 

Component Based Error Estimation. 

For some problem classes, it may be desirable to allow the user to select a particular 
vector component for estimating the error in the solution. A rather short term effort would be 
required to modify the GUI and error estimation algorithms to allow the user to interactively 
select the fifth solution component, for example, on which to estimate the error. 

As indicated above, some of these items are research issues while others are of a develop- 
ment nature. The incorporation of all of these features into Auditor would greatly enhance 
the utility and flexibility of the software package. 
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